From 231a90b9ff564a0328d1050db537902c65e9d19c Mon Sep 17 00:00:00 2001 From: Jeremy Howard Date: Mon, 9 Nov 2026 05:38:55 +0000 Subject: [PATCH] Unify native map execution, conversion and shared array reads --- Cargo.toml | 6 +- DEV.md | 13 +- core/src/agreement.rs | 294 ++++++++++++++++++--------- core/src/array.rs | 3 +- core/src/buffer.rs | 54 ++--- core/src/buffer/output.rs | 121 +++++++++++ core/src/complex/bulk.rs | 25 +++ core/src/convert.rs | 185 ++++++++++++++--- core/src/error.rs | 4 +- core/src/events.rs | 157 +++++--------- core/src/facts.rs | 23 ++- core/src/gather.rs | 20 +- core/src/reduce.rs | 12 +- core/src/search.rs | 15 +- core/src/slice.rs | 2 +- core/src/special.rs | 6 +- core/src/tolerant.rs | 8 +- kernels/src/operation.rs | 73 ++++--- kernels/src/operation/map.rs | 60 +++--- kernels/src/operation/map/complex.rs | 8 +- kernels/src/operation/map/io.rs | 123 +++-------- py/src/explain.rs | 6 +- pyproject.toml | 4 +- src/create.rs | 2 +- src/lib.rs | 26 ++- tests/accuracy.rs | 18 +- tests/linalg.rs | 10 +- tests/numeric.rs | 72 +++++-- tests/slices.rs | 22 +- tests/speed_python.txt | 13 ++ tests/test_speed.py | 32 +-- 31 files changed, 920 insertions(+), 497 deletions(-) create mode 100644 core/src/buffer/output.rs diff --git a/Cargo.toml b/Cargo.toml index 0e8bdb9..b0afc6a 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -7,12 +7,12 @@ version = "0.1.59" edition = "2024" [workspace.dependencies] -rustymath = { path = ".", version = ">=0.1.56, <1", default-features = false } +rustymath = { path = ".", version = ">=0.1.59, <1", default-features = false } rustymath-buffer = { path = "buffer", version = ">=0.1.39, <1" } pyo3 = ">=0.29" pyo3-build-config = ">=0.29" -rustymath-core = { path = "core", version = ">=0.1.56, <1", default-features = false } -rustymath-kernels = { path = "kernels", version = ">=0.1.56, <1" } +rustymath-core = { path = "core", version = ">=0.1.59, <1", default-features = false } +rustymath-kernels = { path = "kernels", version = ">=0.1.59, <1" } rustymath-timing = { path = "timing", version = ">=0.1.56, <1" } half = { version = ">=2.7.1, <3", features = ["num-traits"] } hashbrown = ">=0.15, <1" diff --git a/DEV.md b/DEV.md index 4eab798..fa14e21 100644 --- a/DEV.md +++ b/DEV.md @@ -8,9 +8,9 @@ Rustymath drops precision that numerical code doesn't need, for speed and a smal Keep the documented accuracy bounds that callers use. `nbs/accuracy.qmd` lists every exception in one place. A few ulps more, for real speed or a simpler implementation, become a documented exception there, not a reason to keep slower code. Check a change that alters arithmetic against those contracts, not against bit-identical results. Don't add sign canonicalization, compensated arithmetic or defensive passes without an agreed need. Each shortcut that drops precision is a hypothesis that real uses don't need that precision. When a real use needs digits that a shortcut loses, return that operation to the conservative method, or create an optional parameter. -Rustymath and basedpl have no outside users yet, so treat neither as published. Machinery is shared when it has two users: basedpl and rustymath's Python API. Changing rustymath to give them shared machinery is encouraged. Separate machinery in basedpl and rustymath for the same task is an error, unless Jeremy has asked for them to differ there. +Rustymath is BPL's mathematics and array core, with a useful thin Python interface over the same algorithms. Generic numerical and array algorithms belong here; parsing, language/function semantics, fusion, JIT, batching and Program construction belong to BPL. Choose ownership by the work an algorithm performs, not by whether BPL currently calls a particular convenience wrapper. Consolidate duplicate general algorithms in RM, retaining the better implementation from either project. Rustymath and basedpl have no outside users yet; neither API is treated as published. Discuss public Python removals before making them. -Development and planning run on Syd1 under `/home/ubuntu/aai-ws/{basedpl,rustymath}`. Both projects' `meta/` directories are maintained there and stay outside Git. Start at `../basedpl/meta/open.md` for the shared repair process. Mac copies are snapshots; never sync them over newer Syd1 work. +Read `../basedpl/meta/open.md` for the current working host, authoritative source, builds, plans and evidence. Both projects' `meta/` directories stay outside Git. ## Commands and builds @@ -25,7 +25,7 @@ nbdev-test nbdev-preview ``` -Use `cargo develop` to rebuild and install the Python extension after Rust changes. It builds the `release` profile, which development, tests and timings all use; `fastws.cargo.run_tests(root, ...)` builds it too. Bare `cargo test` uses Cargo's built-in `test` profile, which has `release`'s settings but builds again; pass `--release`. Cargo tests compile missing code but do not install it. An already-running Python process retains its loaded extension. `cargo nextest run --cargo-profile release` runs each test in its own process and every test binary at once. `.config/nextest.toml` reports a test as slow at 5 s, stops it at 30 s, and runs `speed` last with nothing beside it. Its `untimed` and `speed` profiles split the speed test from the rest. basedpl's default test command, `fastws.cargo.run_tests([basedpl, rustymath], profiles=('untimed', 'timed', 'speed'))`, runs each profile in turn, at once in each repository that defines it, so `speed` runs after basedpl's timings, with nothing beside it. +Use `cargo develop` to rebuild and install the Python extension after Rust changes. It builds the `release` profile, which development, tests and timings all use; `fastws.cargo.run_tests(root, ...)` builds it too. Bare `cargo test` uses Cargo's built-in `test` profile, which has `release`'s settings but builds again; pass `--release`. Cargo tests compile missing code but do not install it. An already-running Python process retains its loaded extension. `cargo nextest run --cargo-profile release` runs each test in its own process and every test binary at once. `.config/nextest.toml` reports a test as slow at 5 s, stops it at 30 s, and runs `speed` last with nothing beside it. Its `untimed` and `speed` profiles split the speed test from the rest. The shared default command is `fastws.cargo.run_tests([basedpl, rustymath])`. It uses each repository's default nextest profile, including timings and counts. Explicit profile selection is for focused diagnostics. `cargo nextest run --cargo-profile release` and `pytest -q` include the speed tests, `tests/speed.rs` and `tests/test_speed.py`. Their ordinary runs time each case beside builds and other tests, under a shared hold on the timing lock and on one thread. They also check each case's recorded counts of instructions and allocated bytes for each item, which don't depend on what else the machine runs. A case over its budget is timed again under an exclusive hold where no one else holds the lock, and is reported `unconfirmed`, without failing, where someone does. `RUSTYMATH_COSTS=1` makes either test a calibration run, as "Speed measurements" describes. Use one test binary for a focused check. Run `pytest -q` at the end of a piece of work. Before reporting a failure in untouched code, rebuild with the same profile and rerun it. Builds have no execution timeout. Use the documented managed runners with a 30-second test execution budget. Cargo serializes builds in the same target directory. @@ -85,7 +85,7 @@ A generic function compiles once for each combination of types that callers give ## Code map -The dependency layers are `rustymath-core` → `rustymath-kernels` → the public `rustymath` facade; the facade also uses core directly. Core owns formulas, number facts, storage and geometry together, with no dependency on executable kernels. Kernels owns concrete typed execution. Existing public numerical paths are facade reexports, not duplicate implementations. Files below are relative to `core/src/` unless qualified. +The dependency layers are `rustymath-core` → `rustymath-kernels` → the public `rustymath` facade; the facade also uses core directly. Core owns formulas, number facts, storage and geometry, without depending on executable kernels. Kernels owns concrete typed execution. The facade explicitly exports math/array capabilities, scalar formulas and their public traits/macros; kernel-only helpers stay in the lower crates. BPL owns language evaluation, Program/trace/fusion/JIT policy and diagnostic labels. Files below are relative to `core/src/` unless qualified. | Files | Responsibility | |---|---| @@ -96,13 +96,14 @@ The dependency layers are `rustymath-core` → `rustymath-kernels` → the publi | `slice.rs`, `parallel.rs` | Scalar/slice object interfaces, output storage and elementwise scheduling. | | `reduce.rs`, `stats.rs` | Folds, scans, axis traversal and statistical summaries. | | `facts.rs` | What a planner knows about each function object before it reads an item: the `Range` of its results, from its arguments' ranges by interval rules, whether any result fails, and which argument a constant keeps (`Facts`). Also the one table of which reduction folds each function and which objects take characters, which rustymath, its Python crate and basedpl each expand. | -| `structure.rs` | Structural selection plans. `structure.rs` also holds each structural form's shape rule and the steps that read its result (`take_steps`, `rotate_steps`, `select_steps` and the rest), written once over the `Lengths` trait: `Known` computes on numbers, and `Bounded` on a trace's length variables, recording each decision for later runs to check (`evaluate`). Each concrete plan calls them with `Known`. A recording calls them with `Bounded` for a template, `Steps`, which each run fills at its lengths (`Steps::at`). | +| `structure.rs` | Concrete shape rules and positional plans for indexing, take/drop, rotation, windows, transpose and catenation. Shape functions take `usize` dimensions; `take_steps`, `rotate_steps` and related functions build the `Steps` that read each result. `Walk` and `Slide` recognize shared structural geometry. Symbolic lengths, trace decisions and Program construction belong to BPL. | | `arith.rs`, `integer.rs`, `convert.rs` | Wrapped/checked arithmetic, integer functions and vectorizable conversions. | | `distribution.rs`, `random.rs` | Checked distributions, function objects, sampling, rolls and deals. | | `number_theory.rs` | Primality, factorization, prime enumeration/counting and factorial. | | `field.rs`, `exact.rs`, `linalg/exact.rs`, `polynomial.rs` | Generic field algorithms, rational functions, exact matrices and polynomials. | | `search.rs`, `tolerant.rs`, `gather.rs` | Exact/tolerant search and classification, stable sorting, partial selection, comparisons and copying through conversions. | | `error.rs`, `layout.rs`, `array.rs` | Checked array errors; shapes, axis keys, names and size limits; and the borrowed `ArrayView` over source bindings. | +| `events.rs` | Mathematical routes and one generic capture transport (`Record`, schema-owned rows, `Run`, `Sink`). A producer owns its vocabulary and typed row construction; BPL's tables and language routes live in BPL `src/events.rs`. | | `optimize.rs` | Minimization by five methods behind one `Objective` trait, and solving behind one `Mapping` trait: Brent's root finder and Levenberg–Marquardt for nonlinear equations, and CG, MINRES and GMRES for linear ones. | | `linalg.rs`, `linalg/scalar.rs`, `linalg/factor/` | Packed products, scalar operations and decompositions. | | `signal.rs` | Correlations and convolutions, which slide a kernel along leading axes through `linalg`'s mapped products, and overlap-add. | @@ -250,7 +251,7 @@ The native SVD, Hermitian and general solvers compile everywhere. On macOS the S ## Python boundary -A change that adds a user-facing public Rust function also adds its Python binding, following the existing API's names, argument order and dtypes. Every binding follows the conventions of the Rust core that BasedPL shares. Related bindings agree with each other. A function with NumPy's name needn't behave as NumPy's does. `remainder` and `floor_divide` round through the comparison tolerance where one is set, set functions keep first-occurrence order, orderings put NaN last, and `diagonal` and `trace` take the last two axes, as the matrix functions do. +Expose useful mathematical and array capabilities through a thin Python interface to the shared Rust algorithms. A public Rust helper does not require a matching Python binding. Keep Python argument decoding, defaults, scalar-return policy and result wrappers in `py/src/` or `python/rustymath/`; keep generic mathematics and array algorithms in RM even when a Python convenience is their only direct caller. Related bindings share the core's names, argument order and dtypes. A function with NumPy's name needn't behave as NumPy's does. `remainder` and `floor_divide` round through the comparison tolerance where one is set, set functions keep first-occurrence order, orderings put NaN last, and `diagonal` and `trace` take the last two axes, as the matrix functions do. Python operations are eager. `python/rustymath/create.py` gives fill constructors their names/defaults over one native fill primitive, and names the dtype facts that Rust returns. Sequence generation stays in `create`. Structural functions and bracket indexing in `py/src/structure.rs` delegate selector decoding to `buffer/src/selection.rs`, shared with BasedPL, for the selection and writeback paths. Brackets plan through `structure::plan_parts`, the planner that BasedPL's `⌷` and Python indexing share, with broadcasting, as BasedPL's `.[]` does: index arrays pair by broadcasting, by arithmetic's rule for names. A selector with one axis name selects along the axis with that name (`structure::addressed`). A selection by strides, padding or rotations is a view that shares storage, and a selection by index arrays copies its items (`Array::selected`). `reshape` keeps a view's mapping, which reads the same items in order under any layout. A view whose mapping isn't affine, such as padded windows or a reshape of windows, keeps the items that its first whole read gathers, and later reads use them (`ArrayView::keeps`). An affine view, such as a slice or an affine reshape of one, keeps nothing, and a reader that can't use its strides gathers it again. A small selection that doesn't compose with a view's mapping lists its offsets, and a larger one gathers the view first (`agreement::small_read`, which BasedPL's reads of pending arrays share). `take`, `pad` and `windows` read positions past an end by `structure::Edge`'s modes, which BasedPL's `•edge` names: `take` through `Edge::positions`, and `pad` and padded windows through `Steps::padded`. `structure::Windows` plans windows for `windows` as it does for BasedPL's `↕`. A `fill=` is a `Fill` checked at the array's dtype (`Fill::new`), which the selection's mapping reads (`Mapping::filled`). Assignment detaches shared storage, preserving independent values. `rm.at` combines through the same selection with `Array::combine_at` and its function's `Scatter` kernels, sharing `Array::assign`'s checks. A view retains its whole base until `compact()` copies its logical items, beside any items it has gathered, and even when it is a contiguous subrange. Methods call the same binding functions as package-level operations. `python/rustymath/structure.py` composes conveniences such as trace and repeated differences. Differences run subtraction kernels on two mappings of the input, not gathered copies. `python/rustymath/random.py` subclasses the native generator for distribution conveniences and multivariate-normal sampling. `py/src/random.rs` owns its RNG. Sampling, rolls and deals stay in the shared core. diff --git a/core/src/agreement.rs b/core/src/agreement.rs index b027c65..d16c0d0 100644 --- a/core/src/agreement.rs +++ b/core/src/agreement.rs @@ -20,7 +20,7 @@ pub enum Mapping { Steps { base: usize, steps: Vec }, } -/// What an unmapped position reads. Input plans use the output position; Scalar broadcasts one source item. +/// What an unmapped position reads. Input plans use the output position and are invalid beyond their end. Scalar broadcasts one source item. #[derive(Clone, Debug, PartialEq, Eq)] pub enum Missing { Default, Input(ReadPlan), Scalar(usize), Invalid } @@ -59,23 +59,27 @@ impl bool> RunSink for F { fn run(&mut self, at: #[derive(Clone, Debug, PartialEq, Eq)] pub struct ReadPart { pub end: usize, pub plan: ReadPlan } -/// Data-free positional routes shared by readers, compilation and read-back. +/// Data-free positional routes shared by readers, compilation and read-back. Leaves stay inline; maps and joins share their recursive payloads. #[derive(Clone, Debug, PartialEq, Eq)] pub enum ReadRoute { - Input(usize), Index { base: i64, step: i64 }, Map { input: ReadPlan, map: Mapping, missing: ReadPlan }, Join(Arc<[ReadPart]>), Default, Invalid + Input(usize), Index { base: i64, step: i64 }, Map(Arc), Join(Arc<[ReadPart]>), Default, Invalid } +#[derive(Debug, PartialEq, Eq)] +pub struct ReadMap { pub input: ReadPlan, pub map: Mapping, pub missing: ReadPlan } + #[derive(Clone, Debug, PartialEq, Eq)] -pub struct ReadPlan { route: Arc, items: usize } +pub struct ReadPlan { route: ReadRoute, items: usize } impl ReadPlan { - pub fn input(source: usize, items: usize) -> Self { Self { route: Arc::new(ReadRoute::Input(source)), items } } - pub fn fill(items: usize) -> Self { Self { route: Arc::new(ReadRoute::Default), items } } - pub fn index(base: i64, step: i64, items: usize) -> Self { Self { route: Arc::new(ReadRoute::Index { base, step }), items } } + pub fn input(source: usize, items: usize) -> Self { Self { route: ReadRoute::Input(source), items } } + pub fn fill(items: usize) -> Self { Self { route: ReadRoute::Default, items } } + pub fn index(base: i64, step: i64, items: usize) -> Self { Self { route: ReadRoute::Index { base, step }, items } } fn fallback(missing: Missing, items: usize) -> Self { match missing { - Missing::Input(plan) => plan, Missing::Scalar(source) => Self::input(source, 1).remapped(Mapping::Single, items, Missing::Invalid), - Missing::Default => Self::fill(items), Missing::Invalid => Self { route: Arc::new(ReadRoute::Invalid), items } + Missing::Input(plan) => plan.remapped(Mapping::Linear(1), items, Missing::Invalid), + Missing::Scalar(source) => Self::input(source, 1).remapped(Mapping::Single, items, Missing::Invalid), + Missing::Default => Self::fill(items), Missing::Invalid => Self { route: ReadRoute::Invalid, items } } } pub fn copy(source: usize, source_items: usize, map: Mapping, items: usize, missing: Missing) -> Self { @@ -91,7 +95,7 @@ impl ReadPlan { width = width.checked_add(share).ok_or_else(too_large)?; joined.push(ReadPart { end: width, plan }); } - Ok(Self { route: Arc::new(ReadRoute::Join(joined.into())), items: width.checked_mul(rows).ok_or_else(too_large)? }) + Ok(Self { route: ReadRoute::Join(joined.into()), items: width.checked_mul(rows).ok_or_else(too_large)? }) } pub fn len(&self) -> usize { self.items } pub fn is_empty(&self) -> bool { self.items == 0 } @@ -129,7 +133,8 @@ impl ReadPlan { } true } - ReadRoute::Map { input, map, missing } => { + ReadRoute::Map(read) => { + let ReadMap { input, map, missing } = read.as_ref(); if input.is_empty() { return missing.runs_through(first, stride, len, at, sink); } match stride { 1 => { @@ -156,26 +161,27 @@ impl ReadPlan { } } } - pub fn mapped(input: ReadPlan, map: Mapping, missing: ReadPlan, items: usize) -> Self { + fn mapped(input: ReadPlan, map: Mapping, missing: Missing, items: usize) -> Self { if map.in_order(input.items, &[items]) { return input; } if let Some(plan) = Self::single_map(&input, &map, items) { return plan; } - Self { route: Arc::new(ReadRoute::Map { input, map, missing }), items } + Self { route: ReadRoute::Map(Arc::new(ReadMap { input, map, missing: Self::fallback(missing, items) })), items } } /// `map` over `input`'s own map, composed into one map over `input`'s input. A read of the result then walks the source once, not once for each run of `map`. `None` unless `map` stays within `input`, `input`'s misses read a fill or are invalid, and the composition is structural with no longer tables. fn single_map(input: &ReadPlan, map: &Mapping, items: usize) -> Option { - let ReadRoute::Map { input: inner, map: first, missing } = input.route() else { return None }; + let ReadRoute::Map(read) = input.route() else { return None }; + let ReadMap { input: inner, map: first, missing } = read.as_ref(); let frame = map.shape(items); if map.bounds(&frame)?.end > input.items { return None; } let misses = match missing.route() { ReadRoute::Default => Missing::Default, ReadRoute::Invalid => Missing::Invalid, _ => return None }; let composed = first.compose(map, &input.shape_hint(), &frame)?; - (composed.footprint() <= first.footprint() + map.footprint()).then(|| Self::mapped(inner.clone(), composed, Self::fallback(misses, items), items)) + (composed.footprint() <= first.footprint() + map.footprint()).then(|| Self::mapped(inner.clone(), composed, misses, items)) } - pub fn remapped(&self, map: Mapping, items: usize, missing: Missing) -> Self { Self::mapped(self.clone(), map, Self::fallback(missing, items), items) } + pub fn remapped(&self, map: Mapping, items: usize, missing: Missing) -> Self { Self::mapped(self.clone(), map, missing, items) } /// Calls `visit` with each source that this read reaches and its item count. A plan's routes form a tree, so each place that reads a source visits it once. pub fn try_for_each_source(&self, visit: &mut impl FnMut(usize, usize) -> Result<(), E>) -> Result<(), E> { match self.route() { ReadRoute::Input(source) => visit(*source, self.items), - ReadRoute::Map { input, missing, .. } => { input.try_for_each_source(visit)?; missing.try_for_each_source(visit) } + ReadRoute::Map(read) => { read.input.try_for_each_source(visit)?; read.missing.try_for_each_source(visit) } ReadRoute::Join(parts) => { for part in &**parts { part.plan.try_for_each_source(visit)?; } Ok(()) } ReadRoute::Index { .. } | ReadRoute::Default | ReadRoute::Invalid => Ok(()), } @@ -188,10 +194,10 @@ impl ReadPlan { /// This read with each default read giving `value`. fn defaulted(&self, value: i64) -> Self { match self.route() { - ReadRoute::Map { input, map, missing } => Self::mapped(input.defaulted(value), map.clone(), missing.defaulted(value), self.items), + ReadRoute::Map(read) => Self::mapped(read.input.defaulted(value), read.map.clone(), Missing::Input(read.missing.defaulted(value)), self.items), ReadRoute::Join(parts) => { let parts: Vec<_> = parts.iter().map(|part| ReadPart { end: part.end, plan: part.plan.defaulted(value) }).collect(); - Self { route: Arc::new(ReadRoute::Join(parts.into())), items: self.items } + Self { route: ReadRoute::Join(parts.into()), items: self.items } } ReadRoute::Default => Self::index(value, 0, self.items), ReadRoute::Input(_) | ReadRoute::Index { .. } | ReadRoute::Invalid => self.clone(), @@ -201,13 +207,14 @@ impl ReadPlan { let (source, start, stride) = self.run(&range)?; (stride == 1 || range.len() == 1).then(|| (source, start..start + range.len())) } - pub fn shape_hint(&self) -> Vec { match self.route() { ReadRoute::Map { map, .. } => map.shape(self.items), _ => vec![self.items] } } + pub fn shape_hint(&self) -> Vec { match self.route() { ReadRoute::Map(read) => read.map.shape(self.items), _ => vec![self.items] } } /// The one source that a read over a result of `shape` reaches by an affine map, and that map, where every position stays within the source. `None` for a read that joins sources, fills, or maps otherwise. pub fn positional(&self, shape: &[usize]) -> Option<(usize, Mapping)> { if count(shape) != self.items { return None; } match self.route() { ReadRoute::Input(source) => Some((*source, Mapping::Linear(1))), - ReadRoute::Map { input, map, .. } => { + ReadRoute::Map(read) => { + let ReadMap { input, map, .. } = read.as_ref(); if input.is_empty() || map.bounds(shape)?.end > input.len() { return None; } let inner_shape = input.shape_hint(); let (source, inner) = input.positional(&inner_shape)?; @@ -217,11 +224,21 @@ impl ReadPlan { _ => None, } } - /// One source's affine view. Missing positions and non-affine compositions have no borrowed export. - pub fn affine(&self, shape: &[usize]) -> Option<(usize, usize, Vec)> { - let (source, map) = self.positional(shape)?; - let (base, strides) = map.affine(shape)?; - Some((source, base, strides)) + fn affine_into(&self, shape: &[usize], strides: &mut [isize]) -> Option<(usize, usize, Range)> { + if count(shape) != self.items || shape.len() != strides.len() { return None; } + let (source, map, items) = match self.route() { + ReadRoute::Input(source) => (*source, &Mapping::Linear(1), self.items), + ReadRoute::Map(read) => { + let ReadMap { input, map, .. } = read.as_ref(); + let &ReadRoute::Input(source) = input.route() else { return None }; + (source, map, input.len()) + } + ReadRoute::Join(parts) if parts.len() == 1 => return parts[0].plan.affine_into(shape, strides), + _ => return None, + }; + let base = map.affine_axes(shape, &mut |k, step| strides[k] = step)?; + let bounds = affine_bounds(base, strides, shape)?; + (bounds.end <= items).then_some((source, base, bounds)) } /// The inverse read where this plan visits every item of one source once by affine axes. pub fn inverse(&self, source: usize, shape: &[usize], items: usize) -> Option { @@ -248,17 +265,31 @@ impl ReadPlan { }).collect::>>()?; Some(Mapping::from(steps)) } - /// The items that `range` reads as one run of a source: its source, first position and stride, 0 for one item that every position reads and 1 for consecutive items. `runs` must give the whole range as that one run. + /// A cheap proof that `range` reads one source with stride 0 or 1. This follows route depth and never scans positions; `runs` may enumerate a table just to discover its first span. fn run(&self, range: &Range) -> Option<(usize, usize, isize)> { if range.is_empty() || range.end > self.items { return None; } - let mut found = None; - self.runs(range.clone(), &mut |_, len, run| { - if let Run::Source { source, first, stride } = run && len == range.len() && (len == 1 || matches!(stride, 0 | 1)) { - found = Some((source, first, if len == 1 { 1 } else { stride })); + let n = range.len(); + let (source, first, stride) = match self.route() { + &ReadRoute::Input(source) => (source, range.start, 1), + ReadRoute::Map(read) => { + let ReadMap { input, map, missing } = read.as_ref(); + if input.is_empty() { return missing.run(range); } + let (first, stride) = map.run(range.start, n)?; + let stride = if n == 1 { 1 } else { stride }; + if !matches!(stride, 0 | 1) { return None; } + if first >= input.len() { return missing.run(range); } + let end = first.checked_add(if stride == 0 { 1 } else { n })?; + let (source, at, step) = input.run(&(first..end))?; + (source, at, if stride == 0 { 0 } else { step }) } - false - }); - found + ReadRoute::Join(parts) => { + let (part, at, _, after) = joined_at(parts, range.start)?; + if n > after { return None; } + part.run(&(at..at.checked_add(n)?))? + } + ReadRoute::Index { .. } | ReadRoute::Default | ReadRoute::Invalid => return None + }; + Some((source, first, if n == 1 { 1 } else { stride })) } } @@ -365,46 +396,73 @@ impl<'s> Read<'s> { let end = start.checked_add(if stride == 0 { 1 } else { range.len() })?; Some((binding.borrow(start..end)?, stride as usize)) } - /// `run_view`, where the storage holds `dtype`. - fn storage_run(&self, range: &Range, dtype: crate::buffer::DType) -> Option<(crate::buffer::View<'_>, usize)> { - self.run_view(range).filter(|(view, _)| view.dtype() == dtype) - } /// Borrows a contiguous input range; other reads use scratch at exactly the requested length. pub fn part_view<'a>( &'a self, range: Range, dtype: crate::buffer::DType, checked: bool, scratch: &'a mut crate::buffer::Buffer ) -> Result, Error> { self.range(&range)?; - if let Some((view, stride)) = self.storage_run(&range, dtype) { if stride == 1 || range.len() == 1 { return Ok(view); } } - self.read_into(range, dtype, checked, scratch)?; + let run = self.run_view(&range); + if let Some((view, stride)) = run && view.dtype() == dtype && (stride == 1 || range.len() == 1) { return Ok(view); } + self.write_into(range, dtype, checked, scratch, run)?; Ok(scratch.view()) } + + /// A span at its actual storage dtype and the item count copied. Copied spans require sources of one dtype. + pub fn argument<'a>( + self, range: Range, checked: bool, scratch: &'a mut crate::buffer::Buffer + ) -> Result<(crate::buffer::View<'a>, usize), Error> where 's: 'a { + self.range(&range)?; + if let Some((view, _)) = self.run_view(&range) { return Ok((view, 0)); } + let mut dtype = None; + self.plan.try_for_each_source(&mut |source, _| { + let current = self.sources[source].view.dtype(); + if dtype.is_some_and(|dtype| dtype != current) { return Err(Error::new(Kind::Domain, "read sources have different storage dtypes")); } + dtype = Some(current); + Ok(()) + })?; + let dtype = dtype.ok_or_else(|| Error::new(Kind::Domain, "read has no storage dtype"))?; + let copied = range.len(); + self.write_into(range, dtype, checked, scratch, None)?; + Ok((scratch.view(), copied)) + } + /// A typed span borrowed from storage, or one scalar converted once; `None` means array materialization is needed. + pub fn borrowed(&self, range: Range, checked: bool) -> Result>, Error> { + self.range(&range)?; + if range.is_empty() { return Ok(Some(Arg::Many(&[]))); } + match self.run_view(&range) { Some((view, _)) => borrowed_as(view, checked), None => Ok(None) } + } pub fn part<'a, T: crate::buffer::Cast>(&'a self, range: Range, checked: bool, scratch: &'a mut Vec) -> Result, Error> { + self.part_counted(range, checked, scratch).map(|(items, _)| items) + } + /// The same read and the number of items it materialized in scratch; borrowed runs and scalar conversions count zero. + pub fn part_counted<'a, T: crate::buffer::Cast>( + &'a self, range: Range, checked: bool, scratch: &'a mut Vec + ) -> Result<(Arg<'a, T>, usize), Error> { self.range(&range)?; - if range.is_empty() { return Ok(Arg::Many(&[])); } + if range.is_empty() { return Ok((Arg::Many(&[]), 0)); } if let Some((view, stride)) = self.run_view(&range) { - if view.dtype() != T::DTYPE { - if stride == 0 { return Ok(Arg::One(kept(T::cast(view, 0), checked)?)); } - cast_run(view, stride, range.len(), checked, scratch)?; - return Ok(Arg::Many(scratch)); - } - let items = T::of(view).expect("a typed read selected its dtype"); - return Ok(if stride == 0 { Arg::One(items[0]) } else { Arg::Many(items) }); + if let Some(items) = borrowed_as::(view, checked)? { return Ok((items, 0)); } + cast_run(view, stride, range.len(), checked, scratch)?; + return Ok((Arg::Many(scratch), range.len())); } scratch.clear(); + let copied = range.len(); self.read_onto(range, checked, scratch)?; - Ok(Arg::Many(scratch)) + Ok((Arg::Many(scratch), copied)) } /// Whether `part` reads `range` in place, writing no scratch: one run of a source that holds `T`, or one item that every position reads. pub fn in_place(&self, range: &Range) -> bool { - self.run_view(range).is_some_and(|(view, stride)| stride == 0 || view.dtype() == T::DTYPE) + self.run_view(range).is_some_and(|(view, _)| view.len() == 1 || view.dtype() == T::DTYPE) } /// This read as rows over a result of `shape`, which is rows then items, where every source is bound as the plan reads it, the plan maps one source that holds `T` affinely over `shape`, and each row reads one item or a run of items. `None` for any other read. pub fn rows(&self, shape: [usize; 2]) -> Option> { self.bindings().ok()?; - let (source, base, strides) = self.plan.affine(&shape)?; - let (view, first) = self.sources.get(source)?.affine(base, &strides, &shape)?; - let &[step, inner @ (0 | 1)] = &strides[..] else { return None }; - Some(Rows { items: T::of(view)?, first, step, run: inner == 1 }) + let mut strides = [0; 2]; + let (source, base, bounds) = self.plan.affine_into(&shape, &mut strides)?; + let binding = self.sources.get(source)?; + binding.borrow(bounds)?; + let [step, inner @ (0 | 1)] = strides else { return None }; + Some(Rows { items: T::of(binding.view)?, first: base.checked_sub(binding.start)?, step, run: inner == 1 }) } /// Appends items `range` of the plan onto `output` at `T`, a run at a time (`ReadPlan::runs`): a map over one source walks whole where it can (`walk_onto`), and each other run copies a source's run, gives index values or fills. fn read_onto(&self, range: Range, checked: bool, output: &mut Vec) -> Result<(), Error> { @@ -474,7 +532,8 @@ impl<'s> Read<'s> { } /// The whole source and the mapping of a read that maps one source with no fills, for a reader that walks the mapping itself (`linalg::Mapped`). pub fn mapped_source(&self) -> Option<(crate::buffer::View<'s>, &'s Mapping)> { - let ReadRoute::Map { input, map, .. } = self.plan.route() else { return None }; + let ReadRoute::Map(read) = self.plan.route() else { return None }; + let ReadMap { input, map, .. } = read.as_ref(); if map.fills() { return None; } Some((self.walkable(input, map, self.plan.len())?.view, map)) } @@ -499,10 +558,17 @@ impl<'s> Read<'s> { } pub fn read_into(&self, range: Range, dtype: crate::buffer::DType, checked: bool, output: &mut crate::buffer::Buffer) -> Result<(), Error> { self.range(&range)?; + let run = self.run_view(&range); + self.write_into(range, dtype, checked, output, run) + } + fn write_into( + &self, range: Range, dtype: crate::buffer::DType, checked: bool, + output: &mut crate::buffer::Buffer, run: Option<(crate::buffer::View<'s>, usize)> + ) -> Result<(), Error> { if output.dtype() != dtype { *output = crate::buffer::Buffer::with_capacity(dtype, range.len()); } with_dtype!(dtype, T => { let items = T::vec_mut(output).expect("output has the selected dtype"); - if let Some((view, stride)) = self.run_view(&range) { + if let Some((view, stride)) = run { if view.dtype() != dtype { return cast_run(view, stride, range.len(), checked, items); } let source = T::of(view).expect("a typed read selected its dtype"); items.clear(); @@ -548,6 +614,12 @@ fn joined_at(parts: &[ReadPart], position: usize) -> Option<(&ReadPlan, usize, u let part = parts.get(k)?; Some((&part.plan, row * (part.end - start) + column - start, column - start, part.end - column)) } +fn borrowed_as<'a, T: crate::buffer::Cast>(view: crate::buffer::View<'a>, checked: bool) -> Result>, Error> { + if let Some(items) = T::of(view) { return Ok(Some(Arg::of(items))); } + if view.len() == 1 { return Ok(Some(Arg::One(kept(T::cast(view, 0), checked)?))); } + Ok(None) +} + /// The `len` items of a borrowed run of another dtype, as `T` in `output`: one typed loop over the run ([`Cast::extend_cast`](crate::buffer::Cast::extend_cast)), in place of a routed read for each item. A stride of 0 repeats the run's one item. fn cast_run(view: crate::buffer::View<'_>, stride: usize, len: usize, checked: bool, output: &mut Vec) -> Result<(), Error> { output.clear(); @@ -564,7 +636,7 @@ impl Filled { match missing.route() { ReadRoute::Default => Some(Self::Fill), ReadRoute::Invalid => Some(Self::Never), &ReadRoute::Index { base, step: 0 } => Some(Self::Value(base)), &ReadRoute::Input(source) if missing.len() == 1 => Some(Self::Source(source)), - ReadRoute::Map { input, map: Mapping::Single, .. } => match *input.route() { ReadRoute::Input(source) => Some(Self::Source(source)), _ => None }, + ReadRoute::Map(read) if matches!(&read.map, Mapping::Single) => match *read.input.route() { ReadRoute::Input(source) => Some(Self::Source(source)), _ => None }, _ => None } } @@ -589,10 +661,11 @@ pub enum Period<'a> { } impl Period<'_> { - /// How a reader takes this cycle in runs of up to `run` positions: the period, and how many positions it writes the cycle out over (`tiled_cycle`). A cycle shorter than `LONG_CYCLE`, read in order or by walking entries, is written out over `run` positions, so that a run reads it in place at any phase. A longer walk is gathered once with no repeats. `None` for a longer cycle read in order, which a run reads in place, and for repeats. + /// The period and exact backing length written by `tiled_cycle` for this run size. Long plain cycles return `None`; long walking cycles are stored once. pub fn tiling(&self, run: usize) -> Option<(usize, usize)> { match *self { - Self::Cycles(p) | Self::Walks { len: p, .. } if p < LONG_CYCLE => Some((p, run)), Self::Walks { len, .. } => Some((len, 0)), + Self::Cycles(p) | Self::Walks { len: p, .. } if p > 0 && p < LONG_CYCLE => + Some((p, if p < run { p * (run.div_ceil(p) + 1) } else { p })), Self::Walks { len, .. } => Some((len, len)), Self::Cycles(_) | Self::Repeats(_) => None } } @@ -668,6 +741,22 @@ impl<'a> Walked<'a> { } /// The length of the cycle whose items the mapping's positions read again in order (`Period::Cycles`), where it has one. pub fn cycle(&self) -> Option { *self.cycle.get_or_init(|| match self.mapping.period() { Some(Period::Cycles(p)) => Some(p), _ => None }) } + + /// A walking copy, or one cycle written out to cover runs of up to `cover` items at any phase. + /// `None` only when the mapping has neither a tiling period nor explicit walking entries. + pub fn gathered(&self, input: crate::buffer::View<'_>, cover: usize) -> Option { + with_dtype!(input.dtype(), A => { + let data = A::of(input).expect("the view has its dtype"); + let mut items = Vec::new(); + if tiled_cycle(data, self, cover, &mut items).is_none() { + let Mapping::Steps { steps, .. } = &**self else { return None; }; + let n = Steps::count(steps); + items.reserve(n); + expanded_items(self, data, 0, n, &mut items); + } + Some(A::buffer(items)) + }) + } } impl<'a> From<&'a Mapping> for Walked<'a> { fn from(mapping: &'a Mapping) -> Self { Self::new(Cow::Borrowed(mapping)) } } @@ -684,21 +773,34 @@ impl Mapping { pub fn tabled(&self) -> bool { matches!(self, Self::Steps { steps, .. } if steps.iter().any(|s| matches!(s, Steps::Table(_)))) } /// The writes at the places that a [`Mapping::Steps`] reads, in order, as `From` makes one. `None` for any other mapping. pub fn as_writes(&self) -> Option> { if let Self::Steps { base, steps } = self { Some(Writes::Walk { base: *base, steps }) } else { None } } + /// Affine coefficients for aligned axes, without normalizing or exporting a mapping. + fn affine_axes(&self, shape: &[usize], axis: &mut dyn FnMut(usize, isize)) -> Option { + if shape.contains(&0) { for k in 0..shape.len() { axis(k, 0); } return Some(0); } + if let Self::Steps { base, steps } = self { + if steps.iter().map(Steps::len).ne(shape.iter().copied()) { return None; } + return steps.iter().enumerate().try_fold(*base, |at, (k, s)| { + let (start, step) = if s.len() == 1 { (s.get(0)?, 0) } else { s.affine()? }; + axis(k, step); + Some(at.wrapping_add(start)) + }); + } + if let Self::Linear(period) | Self::Tiled(period) = *self { + if period < count(shape) && period != 1 && !shape.iter().rev().scan(1, |p, &n| { *p *= n; Some(*p) }).any(|p| p == period) { return None; } + } + let (base, mut stride) = (self.get(0)?, count(shape)); + for (k, &len) in shape.iter().enumerate() { + stride /= len; + axis(k, if len > 1 { self.get(stride)? as isize - base as isize } else { 0 }); + } + Some(base) + } + /// This mapping as one entry per axis of `shape`, or `None` where offsets aren't axis-separable. pub fn entries(&self, shape: &[usize]) -> Option<(usize, Vec)> { if shape.contains(&0) { return Some((0, shape.iter().map(|&n| Steps::along(n, 0)).collect())); } if let Self::Steps { base, steps } = self { let (offset, steps) = Steps::per_axis(steps, shape)?; return Some((base + offset, steps)); } - if let Self::Linear(period) | Self::Tiled(period) = *self { - let n = count(shape); - if period < n && period != 1 && !shape.iter().rev().scan(1, |p, &n| { - *p *= n; - Some(*p) - }).any(|p| p == period) { return None; } - } - let base = self.get(0)?; - let steps = shape.iter().zip(strides(shape)) - .map(|(&len, stride)| Some(Steps::Stride { start: 0, len, step: if len > 1 { self.get(stride)? as isize - base as isize } else { 0 } })) - .collect::>()?; + let mut steps = Vec::with_capacity(shape.len()); + let base = self.affine_axes(shape, &mut |k, step| steps.push(Steps::Stride { start: 0, len: shape[k], step }))?; Some((base, steps)) } /// Composes positional reads through the logical axes of shape. @@ -789,8 +891,15 @@ impl Mapping { /// The first offset and stride of items `start..start + n`, where a `Steps` mapping reads one regular run. #[inline] pub fn run(&self, start: usize, n: usize) -> Option<(usize, isize)> { - let Self::Steps { base, steps } = self else { return None }; - Steps::run(*base, steps, start, n) + if n == 0 { return None; } + match *self { + Self::Single => Some((0, 0)), Self::Linear(0) | Self::Tiled(0) => None, + Self::Linear(1) => Some((start, 1)), + Self::Linear(r) if n <= r - start % r => Some((start / r, 0)), + Self::Tiled(p) if n <= p - start % p => Some((start % p, 1)), + Self::Steps { base, ref steps } => Steps::run(base, steps, start, n), + _ => None + } } /// The offset for position `i`, or `None` where it reads a fill. // Per-item loops in basedpl call this for every item. The `Steps` arm calls an out-of-line walk to keep `get` small enough to inline there. `Linear(1)` skips the division. @@ -911,6 +1020,11 @@ impl Mapping { match *self { Self::Single => { f(at, range.len(), (items > 0).then_some((0, 0))); } Self::Linear(0) | Self::Tiled(0) => { f(at, range.len(), None); } + Self::Linear(1) => { + let below = items.saturating_sub(at).min(range.len()); + if below > 0 && !f(at, below, Some((at, 1))) { return; } + if below < range.len() { f(at + below, range.len() - below, None); } + } Self::Linear(r) => { while at < range.end { let (k, end) = (at / r, range.end.min((at / r + 1) * r)); @@ -970,7 +1084,7 @@ impl Agreement { let shape = if left.shape().is_empty() { right.shape() } else { left.shape() }; let len = generated_len(shape)?; let map = |l: &Layout| if len > 1 && l.shape() == shape { Mapping::Linear(1) } else { Mapping::Single }; - return Ok(Self { layout: Layout::from(shape.to_vec()), len, left: map(left), right: map(right) }); + return Ok(Self { layout: Layout::from(shape), len, left: map(left), right: map(right) }); } let [xa, ya] = Self::axis_pairs(left, right); Self::mapped(left, right, &xa, &ya, false) @@ -1106,10 +1220,7 @@ const LINE_BLOCK: usize = 1 << 16; /// The one item that every position in `start..start + n` reads, where `mapping` repeats a single item over all of them. `None` where the positions read more than one item. fn single(mapping: &Mapping, data: &[A], start: usize, n: usize) -> Option { - match *mapping { - Mapping::Single => Some(data[0]), Mapping::Linear(r) if r != 1 && n <= r - start % r => Some(data[start / r]), - Mapping::Steps { .. } => mapping.run(start, n).and_then(|(o, step)| (step == 0).then(|| data[o])), _ => None - } + mapping.run(start, n).and_then(|(o, step)| (step == 0).then(|| data[o])) } /// Items `start..start + n` of an argument that `mapping` lays out, read as `expanded` reads them but always as a slice of all `n` items. A repeated item fills `buf`. @@ -1124,18 +1235,19 @@ pub fn expanded_items<'a, A: Element>(mapping: &Walked, data: &'a [A], start: us /// Items `start..start + n` of an argument that `mapping` lays out, where they need no copy: a single item, a contiguous run, items within one repeat, or items of a cycle: within one period of a mapping that cycles (`Period::Cycles`), and as far as `data` holds a `Tiled` cycle, once or written out past its end (`tiled_cycle`). A `Steps` mapping reads a forward run in place through `run`. pub fn in_place<'a, A: Copy>(mapping: &Walked, data: &'a [A], start: usize, n: usize) -> Option> { - if let Mapping::Steps { .. } = **mapping { - // One check of the walk finds both a repeated item and a forward run. - return match mapping.run(start, n) { - Some((o, 0)) => Some(Arg::One(data[o])), Some((o, 1)) => Some(Arg::Many(&data[o..o + n])), - _ => mapping.cycle().filter(|&p| start % p + n <= p).map(|p| Arg::Many(&data[start % p..start % p + n])) - }; - } - if let Some(x) = single(mapping, data, start, n) { return Some(Arg::One(x)); } - match **mapping { - Mapping::Linear(1) => Some(Arg::Many(&data[start..start + n])), - Mapping::Tiled(m) if start % m + n <= data.len() => Some(Arg::Many(&data[start % m..start % m + n])), _ => None - } + match mapping.run(start, n) { + Some((o, 0)) => return data.get(o).copied().map(Arg::One), + Some((o, 1)) => return data.get(o..o.checked_add(n)?).map(Arg::Many), + _ => {} + } + // A `Tiled` source may already hold several periods (`tiled_cycle`); a cyclic Steps view only borrows within one. + let period = match **mapping { + Mapping::Tiled(p) if p > 0 => p, + Mapping::Steps { .. } => mapping.cycle().filter(|&p| p > 0 && n <= p - start % p)?, + _ => return None + }; + let from = start % period; + data.get(from..from.checked_add(n)?).map(Arg::Many) } /// Items `start..start + n` of an argument that `mapping` lays out over the frame, read `in_place` where they can be. Others copy into `buf`. A repeated item fills each of its periods. A cycled argument copies through `gather::cycled`, and a `Steps` mapping through `walk_plan`, over the walk that its reader keeps (`Walked`). @@ -1147,7 +1259,7 @@ pub fn expanded<'a, A: Element>(mapping: &Walked, data: &'a [A], start: usize, n /// The cycle that `mapping`'s positions read from `data`, written into `buf` as `Period::tiling` says for runs of up to `run` positions, with the `Tiled` mapping that then reads it in place: a cycle read in order copies its first period, and one that walks entries (`Period::Walks`) gathers its period through them. A cycle shorter than its cover is written out past its end, so that a run reads it in place at any phase. `None` where `tiling` gives none, leaving `buf` alone. pub fn tiled_cycle(data: &[A], mapping: &Mapping, run: usize, buf: &mut Vec) -> Option { let period = mapping.period()?; - let (p, cover) = period.tiling(run)?; + let (p, n) = period.tiling(run)?; match period { Period::Walks { base, steps, .. } => { copied(&Walked::from(Mapping::Steps { base, steps: steps.to_vec().into() }), data, 0, p, buf); } _ => { @@ -1156,7 +1268,7 @@ pub fn tiled_cycle(data: &[A], mapping: &Mapping, run: usize, buf: & } } buf.truncate(p); - if p < cover { *buf = buf.repeat(cover.div_ceil(p) + 1); } + if p < n { *buf = buf.repeat(n / p); } Some(Mapping::Tiled(p)) } diff --git a/core/src/array.rs b/core/src/array.rs index a1b0dc2..bbbd968 100644 --- a/core/src/array.rs +++ b/core/src/array.rs @@ -40,7 +40,8 @@ impl<'a> ArrayView<'a> { self.sources.get(source)?.borrow(range) } pub fn affine(&self) -> Option<(View<'a>, usize, Vec)> { - let (source, base, strides) = self.plan.affine(self.layout.shape())?; + let (source, mapping) = self.plan.positional(self.layout.shape())?; + let (base, strides) = mapping.affine(self.layout.shape())?; let (view, base) = self.sources.get(source)?.affine(base, &strides, self.layout.shape())?; Some((view, base, strides)) } diff --git a/core/src/buffer.rs b/core/src/buffer.rs index 8f20c6e..8033e9b 100644 --- a/core/src/buffer.rs +++ b/core/src/buffer.rs @@ -780,6 +780,14 @@ casts!(integers u8, i16, i32, i64); casts!(floats half::f16, f32, f64); casts!(complex Complex32, Complex64); +pub(crate) fn kept_as(view: View<'_>, i: usize) -> Result { + let whole = match view { View::Floats(v) => v.get(i).fract() == 0.0, View::Complexes(v) => v.get(i).re.fract() == 0.0, _ => true }; + let lost = || Error::new(Kind::Domain, format!("a {} item doesn't keep its value as {}", view.dtype().name(), T::DTYPE.name())); + if T::DTYPE.is_integral() && !whole { return Err(lost()); } + T::cast(view, i).map_err(|_| lost()) +} + + impl Buffer { /// Flattened positions whose item differs from its type's fill (zero, false or a space). pub fn nonzero(view: View<'_>) -> Vec { @@ -832,12 +840,7 @@ impl Buffer { pub fn cast(view: View<'_>, dtype: DType) -> Result { Self::cast_mapped(view, &crate::agreement::Mapping::Linear(1), view.len(), dtype) } /// The items of `view` at `dtype`, where every item keeps its value. Items convert as [`Buffer::cast`] converts them, and a float that becomes a Boolean or an integer must be whole. A float rounds to a float of another width. A Domain error where an item doesn't keep its value. pub fn kept(view: View<'_>, dtype: DType) -> Result { - let whole = match view { - View::Floats(v) => (0..v.len()).all(|i| v.get(i).fract() == 0.0), View::Complexes(v) => (0..v.len()).all(|i| v.get(i).re.fract() == 0.0), _ => true - }; - let lost = || Error::new(Kind::Domain, format!("a {} item doesn't keep its value as {}", view.dtype().name(), dtype.name())); - if dtype.is_integral() && !whole { return Err(lost()); } - Self::cast(view, dtype).map_err(|_| lost()) + crate::convert::Conversion { input: view.dtype(), output: dtype, kind: crate::convert::ConversionKind::Keep }.buffer(view, true, crate::tolerant::Tolerance::EXACT) } /// Shared conversion dispatch for compact and mapped logical items. Unselected storage never participates in checked conversion. pub fn cast_mapped(view: View<'_>, mapping: &crate::agreement::Mapping, n: usize, dtype: DType) -> Result { @@ -846,7 +849,8 @@ impl Buffer { let missing = Missing::prototype(view.dtype(), &mut sources); let plan = ReadPlan::copy(0, view.len(), mapping.clone(), n, missing); let mut output = Self::with_capacity(dtype, n); - Read { sources: &sources, plan: &plan }.read_into(0..n, dtype, true, &mut output)?; + crate::convert::Conversion { input: view.dtype(), output: dtype, kind: crate::convert::ConversionKind::Cast } + .map(Read { sources: &sources, plan: &plan }, 0..n, true, &mut output)?; Ok(output) } } @@ -854,37 +858,13 @@ impl Buffer { impl Buffer { /// The items of `view` as integers of `dtype`, as BPL's `ₓ` marks and Python's `astype` convert them: a float takes its floor within `t` and saturates at the range, NaN giving 0 ([`crate::convert::floored`]); an integer wraps to the width; a Boolean gives 0 or 1. Other items and targets convert as [`Buffer::cast`] converts them. pub fn floored(view: View<'_>, dtype: DType, t: impl crate::tolerant::Tolerant) -> Result { - match (view, dtype) { - (View::Floats(values), DType::Integer(width)) => { - Ok(with_width!(width, I => I::buffer(with_floats!(values, |values| values.iter().map(|&x| crate::convert::floored::(x, t)).collect())))) - } - (View::Integers(_) | View::Booleans(_), DType::Integer(_)) => Self::cast_results(view, dtype, false), - _ => Self::cast(view, dtype), - } + let tolerance = crate::tolerant::Tolerance::new(t.value()).expect("a tolerance lies in 0..1"); + crate::convert::Conversion { input: view.dtype(), output: dtype, kind: crate::convert::ConversionKind::Floor }.buffer(view, true, tolerance) } pub fn cast_results(view: View<'_>, dtype: DType, checked: bool) -> Result { - match (view, dtype) { - (View::Floats(values), DType::Integer(width)) => { - Ok(with_width!(width, I => I::buffer(with_floats!(values, |values| values.iter().map(|&x| crate::convert::saturated::(x)).collect())))) - } - (View::Integers(values), DType::Integer(width)) if !checked => { - Ok(with_width!(width, T => with_ints!(values, |values| T::buffer(values.iter().map(|&x| x as T).collect())))) - } - (View::Booleans(values), DType::Integer(width)) => Ok(with_width!(width, T => T::buffer(values.iter().map(|&x| x as T).collect()))), - (_, DType::Character) if view.dtype().is_integral() => { - let items = (0..view.len()).map(|i| { - let value = integer(view, i).ok().and_then(|x| u32::try_from(x).ok()).and_then(char::from_u32); - match value { - Some(value) => Ok(value), None if checked => Err(Error::new(Kind::Domain, "result is not a Unicode character")), None => Ok(char::FILL) - } - }).collect::, _>>()?; - Ok(Self::Characters(items)) - } - _ => with_dtype!(dtype, T => { - let mut items = Vec::new(); - T::extend_cast(view, view.len(), checked, &mut items)?; - Ok(T::buffer(items)) - }), - } + crate::convert::Conversion { input: view.dtype(), output: dtype, kind: crate::convert::ConversionKind::Truncate }.buffer(view, checked, crate::tolerant::Tolerance::EXACT) } } + +mod output; +pub use output::{MapOutput, Output, OutputItem, Scratch}; diff --git a/core/src/buffer/output.rs b/core/src/buffer/output.rs new file mode 100644 index 0000000..9afb20e --- /dev/null +++ b/core/src/buffer/output.rs @@ -0,0 +1,121 @@ +use super::*; +use crate::{Arg, parallel::Slots}; + +pub trait OutputItem: Element { fn slots<'s, 'b>(output: &'s mut MapOutput<'_, 'b>) -> Option<&'s mut Slots<'b, Self>>; } + +macro_rules! outputs { + ([$($v:ident: ($dtype:expr) $t:ty),+]) => { + /// Where a map appends: an owned buffer or the next slots of a caller's result. + pub enum MapOutput<'a, 'b> { Buffer(&'a mut Buffer), $($v(&'a mut Slots<'b, $t>)),+ } + impl<'b> MapOutput<'_, 'b> { + fn dtype(&self) -> DType { match self { Self::Buffer(b) => b.dtype(), $(Self::$v(_) => <$t>::DTYPE),+ } } + pub fn reborrow(&mut self) -> MapOutput<'_, 'b> { match self { Self::Buffer(b) => MapOutput::Buffer(b), $(Self::$v(s) => MapOutput::$v(s)),+ } } + pub fn fill(&mut self, n: usize) { + match self { + Self::Buffer(b) => with_dtype!(b.dtype(), T => { let v = T::vec_mut(b).expect("a buffer has its dtype"); v.resize(v.len() + n, T::FILL); }), + $(Self::$v(s) => { s.pairs(Arg::One(<$t>::FILL), Arg::One(()), n, |x, ()| (x, true)); }),+ + } + } + } + $(impl OutputItem for $t { + fn slots<'s, 'b>(output: &'s mut MapOutput<'_, 'b>) -> Option<&'s mut Slots<'b, Self>> { + match output { MapOutput::$v(s) => Some(s), _ => None } + } + } + impl<'a, 'b> From<&'a mut Slots<'b, $t>> for MapOutput<'a, 'b> { + fn from(slots: &'a mut Slots<'b, $t>) -> Self { MapOutput::$v(slots) } + })+ + }; +} +element_types!((outputs)); + +/// Conversion buffers retained for one execution, by computing dtype and operand position. +#[derive(Default)] +pub struct Scratch { buffers: Vec<(DType, Vec)>, results: Vec, copied: usize } + +impl Scratch { + fn buffers(&mut self, dtype: DType, n: usize) -> &mut [Buffer] { + let at = self.buffers.iter().position(|(stored, _)| *stored == dtype).unwrap_or_else(|| { + self.buffers.push((dtype, Vec::new())); + self.buffers.len() - 1 + }); + let buffers = &mut self.buffers[at].1; + buffers.resize_with(n.max(buffers.len()), || Buffer::with_capacity(dtype, 0)); + &mut buffers[..n] + } + pub fn copied(&self) -> usize { self.copied } + pub fn record_copies(&mut self, items: usize) { self.copied += items; } + fn result(&mut self) -> &mut Vec { + let at = self.results.iter().position(|buffer| buffer.dtype() == T::DTYPE).unwrap_or_else(|| { + self.results.push(Buffer::with_capacity(T::DTYPE, 0)); + self.results.len() - 1 + }); + T::vec_mut(&mut self.results[at]).expect("a result buffer has its computing dtype") + } + pub fn buffer(&mut self, dtype: DType) -> &mut Buffer { &mut self.buffers(dtype, 1)[0] } + pub fn inputs(&mut self) -> [&mut Vec; N] { + let buffers: &mut [Buffer; N] = self.buffers(T::DTYPE, N).try_into().expect("scratch has an entry per operand"); + buffers.each_mut().map(|buffer| T::vec_mut(buffer).expect("scratch keeps its computing dtype")) + } +} + +pub struct Output<'a, 'b> { + pub target: MapOutput<'a, 'b>, pub result_dtype: DType, pub checked: bool, pub scratch: Option<&'a mut Scratch> +} + +impl Output<'_, '_> { + /// The block body is erased before output conversion. + pub fn write( + &mut self, n: usize, + apply: &(dyn Fn(std::ops::Range, &mut Slots<'_, O>, Option<&mut Scratch>) -> bool + Sync) + ) -> Result<(), Error> { + let checked = self.checked; + if let MapOutput::Buffer(buffer) = &mut self.target { + if let Some(values) = O::vec_mut(buffer) { extended(values, n, self.scratch.as_deref_mut(), apply); } else { + with_dtype!(buffer.dtype(), T => { + let values = T::vec_mut(buffer).expect("a buffer has its dtype"); + extended(values, n, self.scratch.as_deref_mut(), &|part, slots, scratch| { + converted::(&mut slots.into(), part, checked, scratch, apply).map_err(crate::error::record).is_ok() + }); + }); + } + return Ok(()); + } + if let Some(slots) = O::slots(&mut self.target) { apply(0..n, slots, self.scratch.as_deref_mut()); return Ok(()); } + converted(&mut self.target, 0..n, checked, self.scratch.as_deref_mut(), apply) + } +} + +fn extended( + values: &mut Vec, n: usize, scratch: Option<&mut Scratch>, + apply: &(dyn Fn(std::ops::Range, &mut Slots<'_, O>, Option<&mut Scratch>) -> bool + Sync) +) { + match scratch { + Some(scratch) => { Slots::append(values, n, |slots| apply(0..n, slots, Some(scratch))); } + None => { crate::parallel::extended(values, n, crate::CHUNK, &|part, slots| apply(part, slots, None)); } + } +} + +fn converted( + target: &mut MapOutput<'_, '_>, range: std::ops::Range, checked: bool, mut work: Option<&mut Scratch>, + apply: &(dyn Fn(std::ops::Range, &mut Slots<'_, O>, Option<&mut Scratch>) -> bool + Sync) +) -> Result<(), Error> { + let mut values = work.as_deref_mut().map_or_else(Vec::new, |work| std::mem::take(work.result::())); + let mut local = Scratch::default(); + let conversion = crate::convert::Conversion { input: O::DTYPE, output: target.dtype(), kind: crate::convert::ConversionKind::Truncate }; + let result = (|| { + for at in (range.start..range.end).step_by(crate::CHUNK) { + let n = crate::CHUNK.min(range.end - at); + values.clear(); + Slots::append(&mut values, n, |slots| apply(at..at + n, slots, work.as_deref_mut())); + if let Err(error) = conversion.append(O::view(&values), checked, crate::tolerant::Tolerance::EXACT, + work.as_deref_mut().unwrap_or(&mut local), target.reborrow()) { + target.fill(range.end - at - n); + return Err(error); + } + } + Ok(()) + })(); + if let Some(work) = work { *work.result::() = values; } + result +} diff --git a/core/src/complex/bulk.rs b/core/src/complex/bulk.rs index a642de3..3e785f0 100644 --- a/core/src/complex/bulk.rs +++ b/core/src/complex/bulk.rs @@ -38,6 +38,31 @@ fn sin_cos(y: &[R], sin: &mut Vec, cos: &mut Vec) { } } +/// Real-input cis through the shared real sin/cos slices, with exact quarter turns. +pub(crate) fn real_cis( + k: &crate::Cis, ops: &NativeOps<{ FUSED }>, [arg]: [Arg<'_, R>; 1], n: usize, out: &mut [MaybeUninit>]) { + let ys = match arg { + Arg::One(y) => { + let value = k.complex(ops, [c(y, R::zero())]).expect("cis has a complex scalar body"); + for slot in &mut out[..n] { slot.write(value); } + return; + } + Arg::Many(ys) => &ys[..n], + }; + let (mut sin, mut cos) = (Vec::new(), Vec::new()); + for (ys, out) in ys.chunks(crate::CHUNK).zip(out[..n].chunks_mut(crate::CHUNK)) { + sin.clear(); + cos.clear(); + sin_cos(ys, &mut sin, &mut cos); + for ((slot, &re), &im) in out.iter_mut().zip(&cos).zip(&sin) { slot.write(c(re, im)); } + if ys.iter().fold(false, |turned, &y| turned | turns(y).1) { + let values = unsafe { slice::written(out) }; + for (value, &y) in values.iter_mut().zip(ys) { if let Some((re, im)) = quarter(y) { *value = c(re, im); } } + } + } +} + + fn split, R: Native + Cast, const N: usize>( k: &K, ops: &NativeOps<{ FUSED }>, args: [Arg<'_, Complex>; N], n: usize, out: &mut [MaybeUninit>], fast: impl FnOnce([(&[R], &[R]); N], &mut Vec, &mut Vec), rare: impl Fn([Complex; N], Complex) -> bool) { diff --git a/core/src/convert.rs b/core/src/convert.rs index ea0184c..91abcb8 100644 --- a/core/src/convert.rs +++ b/core/src/convert.rs @@ -7,11 +7,56 @@ use num_traits::Bounded; pub trait Int: Copy + Sealed + sealed::Ways {} #[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)] -pub enum ConversionKind { Cast, Truncate, Bits, Keep } +pub enum ConversionKind { Cast, Truncate, Bits, Keep, Floor } /// Resolved kinds include the working float width even when the other kind is a word. #[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)] pub struct Conversion { pub input: crate::buffer::DType, pub output: crate::buffer::DType, pub kind: ConversionKind } + +macro_rules! conversion { + ($call:expr, $checked:expr, $tolerance:expr, $write:ident) => {{ + use crate::{buffer::{DType, Element, FloatWidth, Width}, error::{Error, Kind}}; + match ($call.kind, $call.input, $call.output) { + (ConversionKind::Floor, DType::Float(width), DType::Integer(to)) => with_width!(to, I => with_float_width!(width, F => { + crate::with_tolerance!($tolerance, |t| $write!(F, I, |x| Ok(floored(x, t)))) + })), + (ConversionKind::Floor, DType::Integer(from), DType::Integer(to)) => { + with_width!(to, I => with_width!(from, S => $write!(S, I, |x| Ok(x as I)))) + } + (ConversionKind::Floor | ConversionKind::Truncate, DType::Boolean, DType::Integer(to)) => { + with_width!(to, I => $write!(bool, I, |x| Ok(u8::from(x) as I))) + } + (ConversionKind::Truncate, DType::Float(width), DType::Integer(to)) => { + with_width!(to, I => with_float_width!(width, F => $write!(F, I, |x| Ok(saturated(x))))) + } + (ConversionKind::Cast | ConversionKind::Truncate, DType::Integer(from), DType::Integer(to)) => with_width!(to, I => with_width!(from, S => { + $write!(S, I, |x| { + if $checked { I::try_from(x).map_err(|_| Error::new(Kind::Overflow, "integer doesn't fit its converted dtype")) } else { Ok(x as I) } + }) + })), + (ConversionKind::Truncate, input_type, DType::Character) if input_type.is_integral() => $write!(i64, char, |x| { + match u32::try_from(x).ok().and_then(char::from_u32) { + Some(x) => Ok(x), None if $checked => Err(Error::new(Kind::Domain, "result is not a Unicode character")), None => Ok(char::FILL) + } + }), + (ConversionKind::Bits, DType::Float(FloatWidth::F32), DType::Integer(Width::I64)) => { + $write!(f32, i64, |x| Ok(i64::from(x.to_bits()))) + } + (ConversionKind::Bits, DType::Float(FloatWidth::F64), DType::Integer(Width::I64)) => { + $write!(f64, i64, |x| Ok(x.to_bits() as i64)) + } + (ConversionKind::Bits, DType::Integer(_), DType::Float(FloatWidth::F32)) => { + $write!(i64, f32, |x| Ok(f32::from_bits(x as u32))) + } + (ConversionKind::Bits, DType::Integer(_), DType::Float(FloatWidth::F64)) => { + $write!(i64, f64, |x| Ok(f64::from_bits(x as u64))) + } + (ConversionKind::Keep, _, dtype) => with_dtype!(dtype, T => $write!(@keep T)), + (ConversionKind::Cast | ConversionKind::Floor | ConversionKind::Truncate, _, dtype) => with_dtype!(dtype, T => $write!(@cast T)), + _ => $write!(@unsupported), + } + }}; +} impl Conversion { pub fn validity + ?Sized>(self, o: &O, x: O::V) -> Option<(crate::error::Kind, O::B)> { use crate::buffer::{DType, Width}; @@ -45,37 +90,101 @@ impl Conversion { .validity::(&crate::recipe::RangeOps(self.input), crate::domain::Range::of_dtype(self.input)) .is_some_and(|(_, valid)| valid != Some(true))) } - pub fn map( + /// Converts contiguous storage through the shared typed conversion rules, without a read plan. + pub(crate) fn buffer(self, view: crate::buffer::View<'_>, checked: bool, tolerance: crate::tolerant::Tolerance) -> Result { + use crate::buffer::Cast; + macro_rules! write { + ($s:ty, $t:ty, $convert:expr) => { Ok(converted_buffer::<$s, $t>(view, checked, $convert)) }; + (@keep $t:ty) => { Ok(<$t>::buffer((0..view.len()).map(|i| converted_item(crate::buffer::kept_as::<$t>(view, i)).0).collect())) }; + (@cast $t:ty) => {{ + let mut values = Vec::with_capacity(view.len()); + if let Err(error) = <$t>::extend_cast(view, view.len(), checked, &mut values) { crate::error::record(error); values.resize(view.len(), <$t>::FILL); } + Ok(<$t>::buffer(values)) + }}; + (@unsupported) => { Err(crate::error::Error::new(crate::error::Kind::Domain, "unsupported resolved conversion")) }; + } + crate::error::checked(|| conversion!(self, checked, tolerance, write)) + } + pub(crate) fn append<'a, 'b>( + self, view: crate::buffer::View<'_>, checked: bool, tolerance: crate::tolerant::Tolerance, + scratch: &'a mut crate::buffer::Scratch, output: crate::buffer::MapOutput<'a, 'b> + ) -> Result<(), crate::error::Error> { + let sources = [crate::agreement::Binding::full(view)]; + let plan = crate::agreement::ReadPlan::input(0, view.len()); + self.map_in(crate::agreement::Read { sources: &sources, plan: &plan }, 0..view.len(), checked, tolerance, scratch, output) + } + pub(crate) fn map( self, input: crate::agreement::Read<'_>, range: std::ops::Range, checked: bool, output: &mut crate::buffer::Buffer ) -> Result<(), crate::error::Error> { - use crate::buffer::{Buffer, DType, FloatWidth, Floats}; - use crate::error::{Error, Kind}; - if self.kind == ConversionKind::Cast && !matches!((self.input, self.output), (DType::Integer(_), DType::Integer(_))) { - return input.read_into(range, self.output, checked, output); + if output.dtype() != self.output { *output = crate::buffer::Buffer::with_capacity(self.output, range.len()); } else { output.clear(); } + self.map_in(input, range, checked, crate::tolerant::Tolerance::EXACT, &mut crate::buffer::Scratch::default(), crate::buffer::MapOutput::Buffer(output)) + } + + pub fn map_in<'a, 'b>( + self, input: crate::agreement::Read<'_>, range: std::ops::Range, checked: bool, tolerance: crate::tolerant::Tolerance, + scratch: &'a mut crate::buffer::Scratch, target: crate::buffer::MapOutput<'a, 'b> + ) -> Result<(), crate::error::Error> { + let mut output = crate::buffer::Output { target, result_dtype: self.output, checked, scratch: Some(scratch) }; + macro_rules! write { + ($s:ty, $t:ty, $convert:expr) => { written::<$s, $t>(input, range, checked, &mut output, $convert) }; + (@cast $t:ty) => { write!($t, $t, Ok) }; + (@keep $t:ty) => { kept_mapped::<$t>(self.input, input, range, checked, &mut output) }; + (@unsupported) => {{ output.target.fill(range.len()); Err(Error::new(Kind::Domain, "unsupported resolved conversion")) }}; } - let mut scratch = Buffer::with_capacity(self.input, 0); - let values = input.part_view(range, self.input, checked, &mut scratch)?; - *output = match (self.kind, values, self.output) { - (ConversionKind::Cast, _, _) => Buffer::cast_results(values, self.output, checked)?, - (ConversionKind::Keep, _, _) => Buffer::kept(values, self.output)?, - (ConversionKind::Truncate, crate::buffer::View::Floats(_), DType::Integer(_)) => Buffer::cast_results(values, self.output, false)?, - (ConversionKind::Bits, crate::buffer::View::Floats(Floats::F32(values)), DType::Integer(crate::buffer::Width::I64)) => { - Buffer::Integers(crate::buffer::IntBuf::I64(values.iter().map(|x| i64::from(x.to_bits())).collect())) - } - (ConversionKind::Bits, crate::buffer::View::Floats(Floats::F64(values)), DType::Integer(crate::buffer::Width::I64)) => { - Buffer::Integers(crate::buffer::IntBuf::I64(values.iter().map(|x| x.to_bits() as i64).collect())) - } - (ConversionKind::Bits, crate::buffer::View::Integers(values), DType::Float(FloatWidth::F32)) => { - Buffer::Floats(crate::buffer::FloatBuf::F32((0..values.len()).map(|i| f32::from_bits(values.get(i) as u32)).collect())) - } - (ConversionKind::Bits, crate::buffer::View::Integers(values), DType::Float(FloatWidth::F64)) => { - Buffer::Floats(crate::buffer::FloatBuf::F64((0..values.len()).map(|i| f64::from_bits(values.get(i) as u64)).collect())) + crate::error::checked(|| conversion!(self, checked, tolerance, write)) + } +} + +#[inline] +fn converted_item(value: Result) -> (T, bool) { + match value { Ok(value) => (value, true), Err(error) => { crate::error::record(error); (T::FILL, false) } } +} + +fn converted_buffer( + view: crate::buffer::View<'_>, checked: bool, convert: impl Fn(S) -> Result +) -> crate::buffer::Buffer { + let (mut values, mut scratch) = (Vec::with_capacity(view.len()), Vec::new()); + for from in (0..view.len()).step_by(crate::CHUNK) { + let part = view.slice(from..(from + crate::CHUNK).min(view.len())); + let items = if let Some(items) = S::of(part) { items } else { + scratch.clear(); + if let Err(error) = S::extend_cast(part, part.len(), checked, &mut scratch) { + crate::error::record(error); values.resize(values.len() + part.len(), T::FILL); continue; } - (ConversionKind::Truncate, _, _) if self.input.is_integral() && self.output.is_integral() => Buffer::cast_results(values, self.output, false)?, - _ => return Err(Error::new(Kind::Domain, "unsupported resolved conversion")), + &scratch }; - Ok(()) + values.extend(items.iter().map(|&x| converted_item(convert(x)).0)); } + T::buffer(values) +} + +#[inline] +fn written( + input: crate::agreement::Read<'_>, range: std::ops::Range, checked: bool, output: &mut crate::buffer::Output<'_, '_>, + convert: impl Fn(S) -> Result + Sync +) -> Result<(), crate::error::Error> { + output.write::(range.len(), &|part, slots, work| { + let part = part.start + range.start..part.end + range.start; + let in_place = input.in_place::(&part); + let len = if in_place { part.len().max(1) } else { crate::CHUNK }; + let mut local = Vec::new(); + let (mut work, mut copied) = (work, 0); + let scratch = match work.as_deref_mut() { Some(work) if !in_place => { let [scratch] = work.inputs::(); scratch }, _ => &mut local }; + let ok = part.clone().step_by(len).fold(true, |ok, start| { + let range = start..(start + len).min(part.end); + let n = range.len(); + let input = match input.part_counted::(range, checked, scratch) { + Ok((input, n)) => { copied += n; input } + Err(error) => { + crate::error::record(error); + return unsafe { slots.fill(n, |out| { for slot in out { slot.write(T::FILL); } false }) } & ok; + } + }; + slots.pairs(input, crate::Arg::One(()), n, |x, ()| converted_item(convert(x))) & ok + }); + if let Some(work) = work { work.record_copies(copied); } + ok + }) } pub(crate) mod sealed { @@ -226,3 +335,27 @@ ints! { u32: wide, |i| f64::from_low(i.into(), 0.0) as f32, f64::from_low(i.into(), 0.0); i64: wide, |i| wide_float(i) as f32, wide_float(i); } + +fn kept_mapped( + dtype: crate::buffer::DType, input: crate::agreement::Read<'_>, range: std::ops::Range, checked: bool, output: &mut crate::buffer::Output<'_, '_> +) -> Result<(), crate::error::Error> { + output.write::(range.len(), &|part, slots, work| { + let mut local = crate::buffer::Buffer::with_capacity(dtype, 0); + let scratch = match work { Some(work) => work.buffer(dtype), None => &mut local }; + match input.part_view(part.start + range.start..part.end + range.start, dtype, checked, scratch) { + Ok(values) => unsafe { slots.fill(part.len(), |out| { + let mut good = true; + for (i, slot) in out.iter_mut().enumerate() { + slot.write(match crate::buffer::kept_as::(values, i) { + Ok(value) => value, Err(error) => { crate::error::record(error); good = false; T::FILL } + }); + } + good + }) }, + Err(error) => { + crate::error::record(error); + unsafe { slots.fill(part.len(), |out| { for slot in out { slot.write(T::FILL); } false }) } + } + } + }) +} diff --git a/core/src/error.rs b/core/src/error.rs index 9f6cd9c..a5f49d1 100644 --- a/core/src/error.rs +++ b/core/src/error.rs @@ -49,7 +49,7 @@ impl std::fmt::Display for Error { fn fmt(&self, f: &mut std::fmt::Formatter<'_> impl std::error::Error for Error {} -/// Measures and the Interfaces' errors entry: the session's 1 error slot, homed here so kernels set it from every thread that runs the read's work (rustymath's parallel runtime installs the caller's slot on its workers). It keeps the 1st data error and a count of errors recorded. basedpl's `Session` holds one, rustymath's Python API reads its own; its representation is this module's. +/// The first recorded numerical error and its count. The parallel runtime installs the caller's active slot on workers; `checked` supplies one for a standalone operation. #[derive(Clone, Default)] pub struct ErrorSlot(std::sync::Arc); @@ -88,6 +88,8 @@ pub fn record(error: Error) { } /// Checked eager entry for Rust and Python callers. +/// Returns direct failures and, when no slot is active, the first recorded failure. With an +/// active slot, recorded failures remain for its enclosing operation to collect. pub fn checked(run: impl FnOnce() -> Result) -> Result { let (mut run, mut result) = (Some(run), None); let recorded = recorded(&mut || result = run.take().map(|run| run())); diff --git a/core/src/events.rs b/core/src/events.rs index 5f1226f..480adc5 100644 --- a/core/src/events.rs +++ b/core/src/events.rs @@ -1,59 +1,33 @@ -//! Producer events shared by BPL's `•explain` and rustymath's Python `explain`. Capture is process-wide, including workers, with nesting on one thread. Payloads own their diagnostic data and never retain runtime values. -use std::sync::{Arc, Mutex, MutexGuard, PoisonError, atomic::{AtomicUsize, Ordering::Relaxed}}; +//! Diagnostic capture shared by mathematical kernels and their callers. Rows own diagnostic data and never retain runtime values. Capture includes worker threads and permits nesting on one thread. +use std::{borrow::Cow, sync::{Arc, Mutex, MutexGuard, PoisonError, atomic::{AtomicUsize, Ordering::Relaxed}}}; +#[derive(Clone, Copy, Debug, Eq, PartialEq)] +pub struct Route { name: &'static str, taken: bool, declined: bool, itemwise: bool } +impl Route { + pub const fn new(name: &'static str, taken: bool, declined: bool, itemwise: bool) -> Self { Self { name, taken, declined, itemwise } } + pub const fn name(self) -> &'static str { self.name } + pub const fn reports_taken(self) -> bool { self.taken } + pub const fn reports_decline(self) -> bool { self.declined } + pub const fn reports_itemwise(self) -> bool { self.itemwise } +} macro_rules! routes { - ($($id:ident: $name:literal, $taken:literal, $declined:literal, $itemwise:literal;)+) => { - #[derive(Clone, Copy, Debug, Eq, PartialEq)] - pub enum Route { $($id),+ } + ($($id:ident: $name:literal, $declined:literal;)+) => { + #[allow(non_upper_case_globals)] impl Route { pub const ALL: &'static [Self] = &[$(Self::$id),+]; - pub const REQUIRED: &'static [Self] = &[Self::Compiled, Self::FusedFold, Self::Trace]; - pub const fn name(self) -> &'static str { match self { $(Self::$id => $name),+ } } - /// Whether a frontend reports the calls that took this route. - pub const fn reports_taken(self) -> bool { match self { $(Self::$id => $taken),+ } } - /// Whether a frontend reports the calls that declined this route. - pub const fn reports_decline(self) -> bool { match self { $(Self::$id => $declined),+ } } - /// Whether this route handles items one at a time, which basedpl records in `•explain`'s own rows. - pub const fn reports_itemwise(self) -> bool { match self { $(Self::$id => $itemwise),+ } } + $(pub const $id: Self = Self::new($name, true, $declined, false);)+ } }; } routes! { - Trace: "trace", true, false, false; - Fusion: "fusion", true, true, false; - Batching: "batching", true, true, false; - SlidingFold: "sliding fold", true, false, false; - MappedFold: "mapped fold", true, false, false; - AlternatedFold: "alternated fold", true, false, false; - MappedProduct: "mapped product", true, false, false; - VectorKernel: "vector kernel", true, false, false; - MappedKernel: "mapped kernel", true, false, false; - PartialGrade: "partial grade", true, true, false; - PartialSort: "partial sort", true, true, false; - Ranks: "ranks", true, true, false; - KeyClasses: "key classes", true, false, false; - GroupCounts: "group counts", true, false, false; - GroupFolds: "group folds", true, true, false; - Determinant: "determinant", true, true, false; - Permanent: "permanent", true, true, false; - MaskedSelect: "masked select", true, false, false; - Contraction: "contraction", true, true, false; - Program: "program", false, true, false; - FusedFold: "fused fold", true, true, false; - Compiled: "compiled", true, true, false; - InPlace: "in place", true, true, false; - ItemPairs: "item pairs", false, false, true; - ItemLanes: "item lanes", false, false, true; - NumberParts: "number parts", false, false, true; - Accumulation: "accumulation", false, false, true; - NestedLabels: "nested labels", false, false, true; - LabelIds: "label ids", false, false, true; - InverseScan: "inverse scan", false, false, true; - KeyGroups: "key groups", false, false, true; - BatchCells: "batch cells", false, false, true; - BatchedItems: "batched items", false, false, true; - DerivativeItems: "∂ items", false, false, true; - TracedItems: "traced items", false, false, true; + SlidingFold: "sliding fold", false; + MappedFold: "mapped fold", false; + AlternatedFold: "alternated fold", false; + MappedProduct: "mapped product", false; + VectorKernel: "vector kernel", false; + MappedKernel: "mapped kernel", false; + Contraction: "contraction", true; + InPlace: "in place", true; } #[derive(Clone, Copy, Debug, Eq, PartialEq)] @@ -65,79 +39,51 @@ impl EventKind { #[derive(Clone, Debug)] pub struct Record { pub kind: EventKind, pub route: Route, pub reason: &'static str, pub calls: usize, pub items: usize, pub site: Option } -#[derive(Clone, Debug)] -pub enum Event { - Route(Record), Fused { glyphs: String, items: usize, blocks: usize, runs: usize, copied: usize }, - Stored { glyphs: String, items: usize, reason: &'static str, partway: bool }, - Forced { function: String, line: usize, column: usize, items: usize, reason: &'static str, site: Option }, - Program { event: &'static str, steps: usize, listing: String }, Fast { rule: &'static str, function: String }, Fallback { reason: String, items: usize }, - Accumulation { items: usize, size: usize, path: &'static str }, - Mixed { code: String, line: usize, column: usize, target: &'static str, written: &'static str, items: usize }, - Compile { glyphs: String, microseconds: usize }, Itemwise { function: String, line: usize, column: usize, route: Route, items: usize } -} -#[derive(Clone, Copy)] +#[derive(Clone, Copy, Debug, Eq, PartialEq)] pub enum Column { Text, Count, Truth } -pub enum Field<'a> { Text(&'a str), Count(usize), Truth(bool) } +pub type Table = (&'static str, &'static [(&'static str, Column)]); +pub const ROUTES: Table = ("routes", &[("route", Column::Text), ("calls", Column::Count), ("items", Column::Count)]); +pub const DECLINED: Table = ("declined", &[("route", Column::Text), ("reason", Column::Text), ("calls", Column::Count), ("items", Column::Count)]); +pub const STOPPED: Table = ("stopped", &[("reason", Column::Text), ("calls", Column::Count), ("items", Column::Count), ("site", Column::Text)]); + +#[derive(Clone, Debug)] +pub enum Field<'a> { Text(Cow<'a, str>), Count(usize), Truth(bool) } impl Field<'_> { pub fn text(&self) -> &str { match self { Self::Text(t) => t, _ => unreachable!("text column") } } pub fn count(&self) -> usize { match self { Self::Count(n) => *n, _ => unreachable!("count column") } } pub fn truth(&self) -> bool { match self { Self::Truth(t) => *t, _ => unreachable!("truth column") } } + fn borrowed(&self) -> Field<'_> { + match self { Self::Text(t) => Field::Text(t.as_ref().into()), Self::Count(n) => Field::Count(*n), Self::Truth(t) => Field::Truth(*t) } + } } -pub const TABLES: &[(&str, &[(&str, Column)])] = { - use Column::*; - &[ - ("fused", &[("glyphs", Text), ("items", Count), ("blocks", Count), ("runs", Count), ("copied", Count)]), - ("stored", &[("glyphs", Text), ("items", Count), ("reason", Text), ("partway", Truth)]), - ("forced", &[("function", Text), ("line", Count), ("column", Count), ("items", Count), ("reason", Text), ("site", Text)]), - ("programs", &[("event", Text), ("steps", Count), ("listing", Text)]), - ("stopped", &[("reason", Text), ("calls", Count), ("items", Count), ("site", Text)]), ("fast", &[("rule", Text), ("function", Text)]), - ("fallbacks", &[("reason", Text), ("items", Count)]), ("declined", &[("route", Text), ("reason", Text), ("calls", Count), ("items", Count)]), - ("accumulations", &[("items", Count), ("size", Count), ("path", Text)]), - ("mixed", &[("code", Text), ("line", Count), ("column", Count), ("target", Text), ("written", Text), ("items", Count)]), - ("compiles", &[("glyphs", Text), ("microseconds", Count)]), ("routes", &[("route", Text), ("calls", Count), ("items", Count)]), - ("itemwise", &[("function", Text), ("line", Count), ("column", Count), ("route", Text), ("items", Count)]), - ] -}; + +#[derive(Clone, Debug)] +pub enum Event { Route(Record), Row { table: &'static Table, fields: Vec>, label: Option } } impl Event { - pub const fn table(&self) -> usize { + pub fn table(&self) -> &'static Table { match self { - Self::Fused { .. } => 0, Self::Stored { .. } => 1, Self::Forced { .. } => 2, Self::Program { .. } => 3, - Self::Route(r) => match r.kind { EventKind::Stopped => 4, EventKind::Declined => 7, EventKind::Taken => 11 }, Self::Fast { .. } => 5, - Self::Fallback { .. } => 6, Self::Accumulation { .. } => 8, Self::Mixed { .. } => 9, Self::Compile { .. } => 10, Self::Itemwise { .. } => 12 + Self::Route(r) => match r.kind { EventKind::Stopped => &STOPPED, EventKind::Declined => &DECLINED, EventKind::Taken => &ROUTES }, + Self::Row { table, .. } => table, } } - pub const fn name(&self) -> &'static str { TABLES[self.table()].0 } + pub fn name(&self) -> &'static str { self.table().0 } pub fn route(&self) -> Option<&Record> { match self { Self::Route(r) => Some(r), _ => None } } - /// The text that names this event in a route directive, where the event marks a slow path: a stopped trace, a stored or forced pending array, a fallback, a mixed write or an itemwise route. `None` for every other event. + /// The producer's designated text label, if this event has one. pub fn slow_name(&self) -> Option<&str> { match self { - Self::Route(Record { kind: EventKind::Stopped, reason, .. }) | Self::Stored { reason, .. } => Some(*reason), - Self::Forced { function: name, .. } | Self::Fallback { reason: name, .. } | Self::Mixed { code: name, .. } => Some(name.as_str()), - Self::Itemwise { route, .. } => Some(route.name()), _ => None + Self::Route(Record { kind: EventKind::Stopped, reason, .. }) => Some(*reason), + Self::Row { fields, label: Some(i), .. } => Some(fields[*i].text()), _ => None, } } pub fn fields(&self) -> Vec> { use Field::*; match self { Self::Route(r) => match r.kind { - EventKind::Taken => vec![Text(r.route.name()), Count(r.calls), Count(r.items)], - EventKind::Declined => vec![Text(r.route.name()), Text(r.reason), Count(r.calls), Count(r.items)], - EventKind::Stopped => vec![Text(r.reason), Count(r.calls), Count(r.items), Text(r.site.as_deref().unwrap_or_default())], + EventKind::Taken => vec![Text(r.route.name().into()), Count(r.calls), Count(r.items)], + EventKind::Declined => vec![Text(r.route.name().into()), Text(r.reason.into()), Count(r.calls), Count(r.items)], + EventKind::Stopped => vec![Text(r.reason.into()), Count(r.calls), Count(r.items), Text(r.site.as_deref().unwrap_or_default().into())], }, - Self::Fused { glyphs, items, blocks, runs, copied } => vec![Text(glyphs), Count(*items), Count(*blocks), Count(*runs), Count(*copied)], - Self::Stored { glyphs, items, reason, partway } => vec![Text(glyphs), Count(*items), Text(reason), Truth(*partway)], - Self::Forced { function, line, column, items, reason, site } => { - vec![Text(function), Count(*line), Count(*column), Count(*items), Text(reason), Text(site.as_deref().unwrap_or_default())] - } - Self::Program { event, steps, listing } => vec![Text(event), Count(*steps), Text(listing)], - Self::Fast { rule, function } => vec![Text(rule), Text(function)], - Self::Fallback { reason, items } => vec![Text(reason), Count(*items)], - Self::Accumulation { items, size, path } => vec![Count(*items), Count(*size), Text(path)], - Self::Mixed { code, line, column, target, written, items } => { - vec![Text(code), Count(*line), Count(*column), Text(target), Text(written), Count(*items)] - } - Self::Compile { glyphs, microseconds } => vec![Text(glyphs), Count(*microseconds)], - Self::Itemwise { function, line, column, route, items } => vec![Text(function), Count(*line), Count(*column), Text(route.name()), Count(*items)], + Self::Row { fields, .. } => fields.iter().map(Field::borrowed).collect(), } } } @@ -148,9 +94,7 @@ impl Report { pub fn record(&mut self, event: Event) { if let Event::Route(r) = &event { if let Some(existing) = self.records.iter_mut().find_map(|event| match event { - Event::Route(existing) if existing.kind == r.kind && existing.route == r.route && existing.reason == r.reason && existing.site == r.site => { - Some(existing) - } + Event::Route(existing) if existing.kind == r.kind && existing.route == r.route && existing.reason == r.reason && existing.site == r.site => Some(existing), _ => None, }) { existing.calls += r.calls; @@ -164,6 +108,7 @@ impl Report { pub fn routes(&self) -> impl Iterator { self.records.iter().filter_map(Event::route) } } + pub trait Sink: Send + Sync { fn sites(&self) -> bool { false } fn record(&self, event: Event); @@ -212,7 +157,3 @@ pub fn taken(route: Route, items: usize) { emit(|_| Event::Route(Record { kind: pub fn declined(route: Route, reason: &'static str, items: usize) { emit(|_| Event::Route(Record { kind: EventKind::Declined, route, reason, calls: 1, items, site: None })); } -#[inline] -pub fn stopped(reason: &'static str, items: usize, site: impl FnOnce() -> Option) { - emit(|sites| Event::Route(Record { kind: EventKind::Stopped, route: Route::Trace, reason, calls: 1, items, site: if sites { site() } else { None } })); -} diff --git a/core/src/facts.rs b/core/src/facts.rs index a983e26..ad8f561 100644 --- a/core/src/facts.rs +++ b/core/src/facts.rs @@ -188,10 +188,11 @@ impl Op { pub struct NativeSlices { pub(crate) single: Option>, pub(crate) double: Option>, pub(crate) byte: Option>, pub(crate) short: Option>, pub(crate) int: Option>, pub(crate) long: Option>, - pub(crate) complex_single: Option>, pub(crate) complex_double: Option> + pub(crate) complex_single: Option>, pub(crate) complex_double: Option>, + pub(crate) real_complex_single: Option>, pub(crate) real_complex_double: Option> } -/// A slice form takes the item body's run context and writes all `n` results. Real and integer forms return their validity; complex forms leave validation to `complex_valid` at the map boundary. -pub type Slice = fn(&K, &NativeOps<{ FUSED }>, [crate::Arg<'_, T>; N], usize, &mut [std::mem::MaybeUninit]) -> R; +/// A slice form writes all `n` results at `U` from arguments at `T`. The caller reports errors from its verdict `R`: `bool` for real/integer forms, `()` for complex forms. +pub type Slice = fn(&K, &NativeOps<{ FUSED }>, [crate::Arg<'_, T>; N], usize, &mut [std::mem::MaybeUninit]) -> R; impl NativeSlices { pub fn single(&self) -> Option> { self.single } pub fn double(&self) -> Option> { self.double } @@ -201,11 +202,17 @@ impl NativeSlices { pub fn long(&self) -> Option> { self.long } pub fn complex_single(&self) -> Option> { self.complex_single } pub fn complex_double(&self) -> Option> { self.complex_double } + pub fn real_complex_single(&self) -> Option> { self.real_complex_single } + pub fn real_complex_double(&self) -> Option> { self.real_complex_double } - const NONE: Self = Self { single: None, double: None, byte: None, short: None, int: None, long: None, complex_single: None, complex_double: None }; + const NONE: Self = Self { single: None, double: None, byte: None, short: None, int: None, long: None, complex_single: None, complex_double: None, + real_complex_single: None, real_complex_double: None }; fn with_complex(self, single: Slice, double: Slice) -> Self { Self { complex_single: Some(single), complex_double: Some(double), ..self } } + fn with_real_complex(self, single: Slice, double: Slice) -> Self { + Self { real_complex_single: Some(single), real_complex_double: Some(double), ..self } + } } impl NativeSlices where K: Unary + Unary { pub(crate) fn unary() -> Self { @@ -366,6 +373,7 @@ macro_rules! objects { $(native_integer |$nik:pat_param, $nio:tt, $nix:tt| $native_integer:expr;)? $(native_complex |$nfk:pat_param, $nfo:tt, $nfx:tt| $native_complex:expr;)? $(complex_slices $complex_slices:ident;)? + $(real_complex_slices $real_complex_slices:ident;)? $(native_real |$nrk:pat_param, $nro:tt, $nrx:tt| $native_real:expr;)? $(real_complex |$rck:pat_param, $rco:tt, $rcx:tt| $real_complex:expr;)? $(slices $slices:ident;)? @@ -412,6 +420,7 @@ macro_rules! objects { fn native_slices(&self) -> Option> { None $(.or_else(|| Some(NativeSlices::$slices())))? $(.or_else(|| Some(NativeSlices::NONE)).map(|s| s.with_complex(crate::complex::bulk::$complex_slices, crate::complex::bulk::$complex_slices)))? + $(.or_else(|| Some(NativeSlices::NONE)).map(|s| s.with_real_complex(crate::complex::bulk::$real_complex_slices, crate::complex::bulk::$real_complex_slices)))? } } )+}; @@ -552,16 +561,16 @@ objects!(1: crate::tolerant::Floor { formula |k, o, [x]| Some(crate::with_tolerance!(k.0, |t| t.floor_in(o, x))); slope |o, _, _| o.number(F::zero()); integers formula; native_complex |k, o, [x]| Some(crate::tolerant::TolerantComplex::complex_floor(native_tolerance(o.tolerance, k.0), x)); exact |_, [x]| Ok(Some(x.floor())); } crate::tolerant::Ceil { formula |k, o, [x]| Some(crate::with_tolerance!(k.0, |t| t.ceil_in(o, x))); slope |o, _, _| o.number(F::zero()); integers formula; native_complex |k, o, [x]| Some(crate::tolerant::TolerantComplex::complex_ceil(native_tolerance(o.tolerance, k.0), x)); exact |_, [x]| Ok(Some(x.ceil())); } crate::Not { table 1; test |_, o, [x]| Some(o.equal(x, o.number(F::zero()))); integers formula; valid |_, v, o, x, _| truth_valid(v, o, x); } - crate::Cis { output complex; native_complex |_, _, [x]| Some(crate::complex::exp::(num_complex::Complex::new(-x.im, x.re))); complex_slices cis; real_complex |_, o, [x]| (!crate::trig::Kind::Cos.predicate(o, x)).then(|| { let (c, s) = crate::complex::quarter(x).unwrap_or_else(|| (crate::trig::formula(o, x, crate::trig::Kind::Cos), crate::trig::formula(o, x, crate::trig::Kind::Sin))); num_complex::Complex::new(c, s) }); } + crate::Cis { output complex; native_complex |_, _, [x]| Some(crate::complex::exp::(num_complex::Complex::new(-x.im, x.re))); complex_slices cis; real_complex_slices real_cis; real_complex |_, o, [x]| (!crate::trig::Kind::Cos.predicate(o, x)).then(|| { let (c, s) = crate::complex::quarter(x).unwrap_or_else(|| (crate::trig::formula(o, x, crate::trig::Kind::Cos), crate::trig::formula(o, x, crate::trig::Kind::Sin))); num_complex::Complex::new(c, s) }); } crate::LambertW { domain dtype, [x] => crate::lambert::real_domain(dtype, x); native_complex |_, _, [x]| Some(crate::lambert::complex(x)); native_real |_, _, [x]| Some(crate::lambert::lambertw::(x)); checks nan; } crate::Lgamma { slope |o, x, _| o.digamma(x); native_complex |_, _, [x]| Some(crate::complex::lgamma::(x)); native_real |_, _, [x]| Some(crate::special::lgamma::(x)); checks nan; } // ψ's slope is ψ⁽¹⁾. crate::Digamma { slope |o, x, _| o.polygamma(o.number(F::one()), x); native_real |_, _, [x]| Some(crate::special::digamma::(x)); } // 2/√π e^-x², and its negative for erfc. crate::Erf { slope |o, x, _| o.multiply(o.number(F::narrow(std::f64::consts::FRAC_2_SQRT_PI)), o.exp(o.negate(o.multiply(x, x)))); native_real |_, _, [x]| Some(crate::special::erf::(x)); } - crate::Erfc { slope |o, x, _| o.multiply(o.number(F::narrow(-std::f64::consts::FRAC_2_SQRT_PI)), o.exp(o.negate(o.multiply(x, x)))); native_real |_, _, [x]| Some(crate::special::erfc::(x)); } + crate::Erfc { slope |o, x, _| o.multiply(o.number(F::narrow(-std::f64::consts::FRAC_2_SQRT_PI)), o.exp(o.negate(o.multiply(x, x)))); native_real |_, _, [x]| Some(crate::special::erfc::(x)); slices unary; } // √π/2 e^(erfinv x)², from the result. - crate::Erfinv { slope |o, _, r| o.multiply(o.number(F::narrow(std::f64::consts::FRAC_2_SQRT_PI.recip())), o.exp(o.multiply(r, r))); native_real |_, _, [x]| Some(crate::special::erfinv::(x)); checks nan; } + crate::Erfinv { slope |o, _, r| o.multiply(o.number(F::narrow(std::f64::consts::FRAC_2_SQRT_PI.recip())), o.exp(o.multiply(r, r))); native_real |_, _, [x]| Some(crate::special::erfinv::(x)); slices unary; checks nan; } crate::Erfcinv { slope |o, _, r| o.multiply(o.number(F::narrow(-std::f64::consts::FRAC_2_SQRT_PI.recip())), o.exp(o.multiply(r, r))); native_real |_, _, [x]| Some(crate::special::erfcinv::(x)); slices unary; checks nan; } ); objects!(2: diff --git a/core/src/gather.rs b/core/src/gather.rs index f7a8e08..108e25c 100644 --- a/core/src/gather.rs +++ b/core/src/gather.rs @@ -762,7 +762,20 @@ fn stepping(values: &[T], strides: &mut [Steps], position: impl Fn(T) - Some(first) } -impl PartialEq for Table { fn eq(&self, other: &Self) -> bool { self.len == other.len && (0..self.len).all(|j| self.get(j) == other.get(j)) } } +impl PartialEq for Table { + fn eq(&self, other: &Self) -> bool { + if self.len != other.len { return false; } + // Cloned geometry compares without reading its index data; distinct recipes still compare their offsets. + let shared = self.start == other.start && match (&self.all, &other.all) { + (Positions::Listed(a), Positions::Listed(b)) => Arc::ptr_eq(a, b), + (Positions::Indices { values: a, size: sa, scale: ka, .. }, Positions::Indices { values: b, size: sb, scale: kb, .. }) => { + sa == sb && ka == kb && Arc::ptr_eq(a, b) + } + _ => false, + }; + shared || (0..self.len).all(|j| self.get(j) == other.get(j)) + } +} impl Eq for Table {} /// Position `j` of the index values `view` along an axis of `size` positions. @@ -991,13 +1004,14 @@ impl Steps { } Some(offset) } - /// An offset that no item of `steps` passes from `base`: `base` plus each entry's furthest offset, where a fill reaches none and a padded entry stops at the end of its axis. Only a table's entry reads its offsets to find it. `None` where an entry has no items, so the walk reads nothing. + /// An offset that no item of `steps` passes from `base`: `base` plus each entry's furthest offset, where a fill reaches none and a padded entry stops at the end of its axis. Indexed tables use their axis bounds; listed tables read their offsets. `None` where an entry has no items, so the walk reads nothing. pub fn furthest(base: usize, steps: &[impl std::borrow::Borrow]) -> Option { steps.iter().try_fold(base, |at, s| { let furthest = match *s.borrow() { ref s if s.len() == 0 => return None, Self::Stride { start, len, step } => start.wrapping_add_signed(step.max(0) * (len - 1) as isize), Self::Padded { size, stride, .. } => size.saturating_sub(1) * stride, Self::Rotated { len, stride, .. } => (len - 1) * stride, - Self::Table(ref t) => (0..t.len()).filter_map(|j| t.get(j)).max().unwrap_or(0) + Self::Table(ref t) => t.span().map_or_else( + || (0..t.len()).filter_map(|j| t.get(j)).max().unwrap_or(0), |(scale, size)| size.saturating_sub(1) * scale) }; Some(at + furthest) }) diff --git a/core/src/reduce.rs b/core/src/reduce.rs index 6af71f9..fe26e51 100644 --- a/core/src/reduce.rs +++ b/core/src/reduce.rs @@ -523,7 +523,7 @@ mod lanes { for (o, x) in outs.iter_mut().zip(chunks) { let scanned = |k: usize| { // SAFETY: the load reads four items of one chunk. - let v = unsafe { load(x.as_ptr().add(4 * k)) }; + let v = load(x.as_ptr().add(4 * k)); let v = add(v, up1(fill, v)); add(v, up2(fill, v)) }; @@ -533,12 +533,10 @@ mod lanes { let s012 = add(s01, last(p2)); let o = o.as_mut_ptr().cast::(); // SAFETY: each store writes four slots of one chunk. - unsafe { - store(o, add(c, p0)); - store(o.add(4), add(add(c, s0), p1)); - store(o.add(8), add(add(c, s01), p2)); - store(o.add(12), add(add(c, s012), p3)); - } + store(o, add(c, p0)); + store(o.add(4), add(add(c, s0), p1)); + store(o.add(8), add(add(c, s01), p2)); + store(o.add(12), add(add(c, s012), p3)); c = add(c, add(s012, last(p3))); } let mut total = first(c); diff --git a/core/src/search.rs b/core/src/search.rs index f9a15ce..774c1f2 100644 --- a/core/src/search.rs +++ b/core/src/search.rs @@ -1,7 +1,7 @@ //! Search and sort over slices of keys. A key is an item read as a whole number (`Key`), or an unsigned word whose order is the order you want (`RadixKey`). `Float::order_key` gives each float such a word, as wide as the float. Positions take 32 bits inside these functions, so a slice holds at most `u32::MAX` items. `tolerant` builds its searches over floats on these functions. //! -//! `index_of`, `classes` and `firsts` index the first position of each key. `index_of` passes over `x` once for each item of `y` when `y` has at most `SCANNED` items. The index is a table when the keys span at most twice the slice's length, or 1024, and a hash map otherwise. `numbered` numbers the classes in order of first position, with each one's size. `grouped_keys` numbers them the same way from a table of keys, with each one's positions, `numbered_keys` with each item's number, and `class_counts` gives only their sizes. `counts` counts the items with each key in such a table. `sort_rows` sorts keys that span such a table by counting them. It sorts other keys with an LSD radix sort that skips any byte that is the same in every key. `view_sort_along`, `view_partition_along` and `view_quantiles_along` work on each item's `SortKey`, which reads back to the item exactly. A sort counts, sorts a float lane through a bitonic network on vectors, or sorts the keys by comparison or by the same radix sort. With several threads, the radix sort splits each pass of more than `RADIX_PIECE` keys across them, and `index_of` splits its lookups. -use crate::hashing::{HashMap, HashMapExt}; +//! `index_of` and `classes` index the first position of each key; `firsts` marks first occurrences directly. `index_of` passes over `x` once for each item of `y` when `y` has at most `SCANNED` items. The index is a table when the keys span at most twice the slice's length, or 1024, and a hash table otherwise. `numbered` numbers the classes in order of first position, with each one's size. `grouped_keys` numbers them the same way from a table of keys, with each one's positions, `numbered_keys` with each item's number, and `class_counts` gives only their sizes. `counts` counts the items with each key in such a table. `sort_rows` sorts keys that span such a table by counting them. It sorts other keys with an LSD radix sort that skips any byte that is the same in every key. `view_sort_along`, `view_partition_along` and `view_quantiles_along` work on each item's `SortKey`, which reads back to the item exactly. A sort counts, sorts a float lane through a bitonic network on vectors, or sorts the keys by comparison or by the same radix sort. With several threads, the radix sort splits each pass of more than `RADIX_PIECE` keys across them, and `index_of` splits its lookups. +use crate::hashing::{HashMap, HashMapExt, HashSet, HashSetExt}; use crate::{ Float, Native, agreement::Mapping, array::ArrayView, buffer::{Buffer, DType, Element, Floats, IntBuf, View, Width}, error::{Error, Kind}, interrupt::check, layout::{Layout, count}, tolerant::{self, Tolerance} @@ -166,10 +166,10 @@ pub fn classes(x: &[K]) -> Vec { x.iter().map(move |&n| first[n.key().wrapping_sub(min) as usize]).collect() } -/// Whether each item of `x` is the first with its key. The pass that writes each key's first position marks it as it goes, with no second pass over `x`. Wider keys mark the items that are their own class. +/// Whether each item is the first with its key. Narrow key ranges mark first positions in a table; wider keys use a hash set. #[inline(never)] pub fn firsts(x: &[K]) -> Vec { - let Some((min, mut table)) = table(x, u32::MAX) else { return first_mask(&first_classes(x, |n| n)) }; + let Some((min, mut table)) = table(x, u32::MAX) else { return first_mask(x, |n| n) }; let mut mask = vec![false; x.len()]; write_first(x, min, &mut table, None, |p, new| mask[p] |= new); mask @@ -590,8 +590,11 @@ impl Ord for Exact<'_> { impl PartialOrd for Exact<'_> { fn partial_cmp(&self, other: &Self) -> Option { Some(self.cmp(other)) } } -/// Whether each position of `classes` is its own class, which makes its item the first with its key. -pub(crate) fn first_mask(classes: &[u32]) -> Vec { classes.iter().enumerate().map(|(i, &c)| c as usize == i).collect() } +/// Whether each item is the first with its key, from one hash insert per item. +pub(crate) fn first_mask(x: &[T], key: impl Fn(T) -> K) -> Vec { + let mut seen = HashSet::with_capacity(x.len()); + x.iter().map(|&v| seen.insert(key(v))).collect() +} /// The positions of `groups` in order of group, keeping their order within each group, with the number in each of the `n` groups. #[inline(never)] diff --git a/core/src/slice.rs b/core/src/slice.rs index 9c05c63..8f8e1cd 100644 --- a/core/src/slice.rs +++ b/core/src/slice.rs @@ -40,7 +40,7 @@ impl<'a, T: Copy> Arg<'a, T> { /// Items `start..start + n` of this argument. A `One` argument stays as it is. #[inline(always)] - pub(crate) fn part(self, start: usize, n: usize) -> Self { match self { Arg::One(v) => Arg::One(v), Arg::Many(s) => Arg::Many(&s[start..start + n]) } } + pub fn part(self, start: usize, n: usize) -> Self { match self { Arg::One(v) => Arg::One(v), Arg::Many(s) => Arg::Many(&s[start..start + n]) } } } /// A type that function objects take or give: a float, a complex number, an integer, a Boolean, a character or an ordering. Pieces on several threads read and write items. diff --git a/core/src/special.rs b/core/src/special.rs index 2d1f14b..7b2d571 100644 --- a/core/src/special.rs +++ b/core/src/special.rs @@ -662,9 +662,11 @@ fn erfc_blocks(x: &[R], out: &mut [MaybeUninit] true }) } -/// [`erfinv`](crate::erfinv) as a function object. +/// [`erfinv`](crate::erfinv) as a function object, with the slice form `many`, which runs `regions` on a slice of f64s. #[derive(Clone, Copy, Debug)] pub struct Erfinv; +forms!(Erfinv at real: |_, [p]| erfinv::<_, { crate::FUSED }>(p), + many |_, [p], n, out| inverses(p, n, out, #[inline(always)] |p| erfinv::<_, { crate::FUSED }>(p), #[inline(always)] |p| (p.abs(), 1.0 - p.abs(), p))); /// [`erfcinv`](crate::erfcinv) as a function object, with the slice form `many`, which runs `regions` on a slice of f64s. #[derive(Clone, Copy, Debug)] @@ -672,7 +674,7 @@ pub struct Erfcinv; forms!(Erfcinv at real: |_, [q]| erfcinv::<_, { crate::FUSED }>(q), many |_, [q], n, out| inverses(q, n, out, #[inline(always)] |q| erfcinv::<_, { crate::FUSED }>(q), #[inline(always)] |q| ((q - 1.0).abs(), q.min(2.0 - q), 1.0 - q))); -/// The slice form of [`Erfcinv`]: `regions` on a slice of f64s, and `f` on each item otherwise. `parts` gives an f64 item's `maddock` arguments and the value whose sign its result takes. +/// The slice form of [`Erfinv`] and [`Erfcinv`]: `regions` on a slice of f64s, and `f` on each item otherwise. `parts` gives an f64 item's `maddock` arguments and the value whose sign its result takes. fn inverses(x: Arg<'_, R>, n: usize, out: &mut [MaybeUninit], f: impl Fn(R) -> R, parts: impl Fn(f64) -> (f64, f64, f64)) -> bool { match x { Arg::Many(x) if R::WIDE => regions::(&x[..n], out, f, parts), x => each!([x], n, out, |[x]| f(x)) } } diff --git a/core/src/tolerant.rs b/core/src/tolerant.rs index 404cc78..3acb5f8 100644 --- a/core/src/tolerant.rs +++ b/core/src/tolerant.rs @@ -117,6 +117,12 @@ pub trait Tolerant: Copy + Send + Sync { let whole = x.fract() == 0.0 && y.fract() == 0.0 && x * y < 1e14 && self.value() <= Tolerance::DEFAULT.value(); if let Some(g) = whole.then(|| crate::integer::gcd(int(x)?, int(y)?)).flatten() { return float(g); } let (t, (nx, dy)) = (self.value(), ratio_parts(x, y)); + if t == 0.0 { + // Odd significands share no powers of two; the smaller exponent supplies their common unit. + let g = crate::integer::gcd(nx.0 as i64, dy.0 as i64).unwrap() as f64; + let unit = if nx.1 == 0 { x / nx.0 as f64 } else { y / dy.0 as f64 }; + return F::narrow(unit * g); + } if let Some(g) = fixed_gcd(x, y, nx, dy, t) { return F::narrow(g); } // The same walk on big integers, where the test `|p·d − q·n| ≤ t·n·q` is exact for `t = m·2^e`. let (n, d) = (BigInt::from(nx.0) << nx.1, BigInt::from(dy.0) << dy.1); @@ -678,7 +684,7 @@ pub fn classes(t: impl Tolerant, x: &[F]) -> Result, Error> { #[inline(never)] pub fn firsts(t: impl Tolerant, x: &[F]) -> Result, Error> { if let Some((keys, _)) = wholes(t, x, whole_limit(t)) { return Ok(search::firsts(&keys)); } - if t.value() == 0.0 { return Ok(search::first_mask(&search::first_classes(x, F::order_key))); } + if t.value() == 0.0 { return Ok(search::first_mask(x, F::order_key)); } let mut firsts = vec![false; x.len()]; checked_runs(t, x, |run| firsts[run[0]] = true)?; Ok(firsts) diff --git a/kernels/src/operation.rs b/kernels/src/operation.rs index d398951..7fd0312 100644 --- a/kernels/src/operation.rs +++ b/kernels/src/operation.rs @@ -1,7 +1,7 @@ //! Resolved numerical operations over shared [`Read`](crate::agreement::Read) plans. [`Kernel::try_resolve`] fixes the computing dtype, result dtype and character policy from an object's [`Formula`](crate::facts::Formula). The same definition supplies native execution and graph formulas. //! //! [`DType::computes_in`] promotes unsupported Boolean and integer inputs to i64, then to the NumPy float width. Unsupported character arithmetic is a domain error. -use crate::buffer::{Buffer, DType, Element, FloatWidth, Source, View, Width}; +use crate::buffer::{Buffer, DType, Element, FloatWidth, MapOutput, Output, Source, View, Width}; use crate::domain::Range; use crate::error::{Error, Kind}; use crate::gather::Writes; @@ -63,6 +63,28 @@ impl std::fmt::Debug for Kernel { } } +/// A resolved elementwise call. Kind selection is complete before this is executed. +#[derive(Clone, Debug)] +pub enum Map { Kernel(Kernel, Settings), Convert(crate::convert::Conversion, Settings, bool) } +impl Map { + pub fn arity(&self) -> usize { match self { Self::Kernel(kernel, _) => kernel.arity(), Self::Convert(..) => 1 } } + /// Appends every result, including initialized fallback items on a numerical error. + /// Errors are execution failures, not a pre-write refusal to handle the call. + /// Both variants follow [`crate::error::checked`]: an enclosing checked scope owns recorded + /// numerical failures. Under an active error slot, `Ok` alone does not establish success. + /// Discard private outputs when that scope fails; neither `Err` nor a recorded failure is a route refusal. + pub fn run<'a, 'b>( + &self, inputs: &[crate::agreement::Read<'_>], range: std::ops::Range, + scratch: &'a mut crate::buffer::Scratch, mut output: MapOutput<'a, 'b> + ) -> Result<(), Error> { + if inputs.len() != self.arity() { output.fill(range.len()); return Err(Error::new(Kind::Length, "map arity does not match its operands")); } + match self { + Self::Kernel(kernel, settings) => kernel.map_in(inputs, range, *settings, scratch, output), + Self::Convert(conversion, settings, checked) => conversion.map_in(inputs[0], range, *checked, settings.tolerance, scratch, output) + } + } +} + trait Application: Send + Sync { fn exact(&self, args: &[num_rational::BigRational]) -> Result, Error>; fn has_exact(&self) -> bool; @@ -361,7 +383,8 @@ impl Kernel { // f64 endpoints cannot prove i64 wrapping bounds. Narrower integers are exact here until their result exceeds their type. if self.variant.dtype.is_integral() && !self.variant.fails_on_overflow() && self.result.is_integral() && (self.result == DType::Integer(Width::I64) || !range.within(bounds.lo, bounds.hi)) { bounds } - else { range } + // Checked arithmetic can rerun wider; only a Boolean result is bounded independently of that attempt. + else if self.result == DType::Boolean { range.meet(bounds) } else { range } } pub fn identity(&self) -> Option { self.definition.identity() } /// The argument that a dyadic call returns unchanged where `constants` gives the arguments that are one number, 0 for the left and 1 for the right, from the formula's identities: `x+0` and `0+x` keep `x`, `x-0` keeps `x`, and `0-x` keeps nothing. @@ -428,40 +451,36 @@ impl Kernel { } pub fn map(&self, inputs: &[crate::agreement::Read<'_>], range: std::ops::Range, settings: Settings, output: &mut Buffer) -> Result<(), Error> { output.clear(); - self.mapping(Inputs::Mapped(inputs), range, settings, MapOutput::Buffer(output)) + self.mapping(Inputs(inputs), range, settings, MapOutput::Buffer(output), None) } - /// Appends n results from ready inputs into output. A one-item view broadcasts; every other view contains n items. - pub fn map_ready(&self, inputs: &[View<'_>], n: usize, settings: Settings, mut output: MapOutput<'_, '_>) -> Result<(), Error> { - if inputs.iter().any(|v| v.len() != 1 && v.len() != n) { output.fill(n); return Err(Error::new(Kind::Length, "ready kernel input length changed")); } - self.mapping(Inputs::Ready(inputs), 0..n, settings, output) + /// Appends a mapped range using the caller's execution-local conversion buffers. + pub fn map_in<'a, 'b>( + &self, inputs: &[crate::agreement::Read<'_>], range: std::ops::Range, settings: Settings, + scratch: &'a mut crate::buffer::Scratch, output: MapOutput<'a, 'b> + ) -> Result<(), Error> { + self.mapping(Inputs(inputs), range, settings, output, Some(scratch)) } - fn mapping(&self, inputs: Inputs<'_, '_>, range: std::ops::Range, settings: Settings, mut target: MapOutput<'_, '_>) -> Result<(), Error> { + fn mapping<'a, 'b>( + &self, inputs: Inputs<'_, '_>, range: std::ops::Range, settings: Settings, mut target: MapOutput<'a, 'b>, + scratch: Option<&'a mut crate::buffer::Scratch> + ) -> Result<(), Error> { use crate::characters::{ByKinds, Policy}; if inputs.len() != self.arity() { target.fill(range.len()); return Err(Error::new(Kind::Length, "kernel operand count changed")); } if range.is_empty() { return Ok(()); } - let mut output = Output { target, dtype: self.result, checked: self.variant.fails_on_overflow() }; + let mut output = Output { target, result_dtype: self.result, checked: self.variant.fails_on_overflow(), scratch }; crate::error::checked(|| match self.policy { Some(Policy::Kinds(ByKinds::Truth(value))) => map::constant(value, range.len(), &mut output), Some(Policy::Kinds(side)) => map::copied(inputs.selected(usize::from(side == ByKinds::Right)), range, self.variant.checked, &mut output), _ => { - let apply = |inputs: Inputs<'_, '_>, output: &mut Output<'_, '_>| { - self.definition.map(self.variant, inputs, range.clone(), settings, self.repair, self.may_raise(), output) - }; - match inputs { - Inputs::Mapped(reads) => { - self.projected(reads, Some(range.clone()), Seed::Identity, settings, |reads, _| apply(Inputs::Mapped(reads), &mut output)) - } - Inputs::Ready(_) if self.projection(settings) == Projection::AsIs => apply(inputs, &mut output), - Inputs::Ready(views) => { - let projected = match views.iter().map(|&v| projected_view(v, self.projection(settings))).collect::, _>>() { - Ok(projected) => projected, Err(error) => { output.target.fill(range.len()); return Err(error); } - }; - let views: Vec<_> = projected.iter().map(|v| v.view()).collect(); - apply(Inputs::Ready(&views), &mut output) - } - } + let mut entered = false; + let result = self.projected(inputs.0, Some(range.clone()), Seed::Identity, settings, |reads, _| { + entered = true; + self.definition.map(self.variant, Inputs(reads), range.clone(), settings, self.repair, self.may_raise(), &mut output) + }); + if !entered { output.target.fill(range.len()); } + result } }) } @@ -842,5 +861,5 @@ mod binary; mod fold; mod map; mod product; -pub use {at::*, map::MapOutput, product::*}; -use {binary::*, fold::*, map::{Inputs, Output, map_formula}}; +pub use {at::*, product::*}; +use {binary::*, fold::*, map::{Inputs, map_formula}}; diff --git a/kernels/src/operation/map.rs b/kernels/src/operation/map.rs index 794329c..bc15fbf 100644 --- a/kernels/src/operation/map.rs +++ b/kernels/src/operation/map.rs @@ -6,20 +6,15 @@ mod integer; mod real; use {complex::*, integer::*, real::map_real}; mod io; -pub use io::MapOutput; -use io::OutputItem; -pub(super) use io::{Inputs, Output}; +use crate::buffer::{Output, OutputItem}; +pub(super) use io::Inputs; pub(super) fn map_formula, const N: usize>( object: &K, variant: crate::domain::Variant, inputs: Inputs<'_, '_>, range: std::ops::Range, settings: Settings, repair: bool, judged: bool, output: &mut Output<'_, '_>) -> Result<(), Error> { if inputs.len() != N { return Err(Error::new(Kind::Length, "kernel operand count changed")); } let ops = settings.ops(); crate::events::taken(crate::events::Route::VectorKernel, range.len()); - match inputs { - Inputs::Mapped(reads) => report_maps(reads.iter().map(|read| read.plan), range.clone()), - Inputs::Ready(views) if range.len() > 1 && views.iter().any(|v| v.len() == 1) => crate::events::taken(crate::events::Route::MappedKernel, range.len()), - _ => {} - } + report_maps(inputs.0.iter().map(|read| read.plan), range.clone()); // A kernel whose arguments' types can't make an error judges no item: `judge` is `variant` unchecked. Reads and `IntegerOps` keep `variant`, so an argument that doesn't convert, or an overflow that `IntegerOps::result` reports, still fails. let judge = crate::domain::Variant { checked: variant.checked && judged, ..variant }; let slices = object.native_slices(); @@ -50,10 +45,10 @@ pub(super) fn map_formula, const N: usize>( }, DType::Integer(width) => with_width!(width, I => map_integer::(object, variant, inputs, range, judge, output)), DType::Complex(crate::buffer::ComplexWidth::F32) if K::REAL_COMPLEX && inputs.real() => { - map_real_complex::(object, variant, inputs, range, &ops, output) + map_real_complex::(object, variant, inputs, range, &ops, slices.and_then(|s| s.real_complex_single()), output) } DType::Complex(crate::buffer::ComplexWidth::F64) if K::REAL_COMPLEX && inputs.real() => { - map_real_complex::(object, variant, inputs, range, &ops, output) + map_real_complex::(object, variant, inputs, range, &ops, slices.and_then(|s| s.real_complex_double()), output) } DType::Complex(crate::buffer::ComplexWidth::F32) => { map_complex::(object, variant, inputs, range, &ops, slices.and_then(|s| s.complex_single()), judge, output) @@ -112,26 +107,37 @@ fn map_blocks( apply: &(dyn Fn([Arg<'_, F>; N], usize, &mut [std::mem::MaybeUninit]) -> Result<(), Kind> + Sync)) -> Result<(), Error> { let recorded = |result: Result<(), Kind>| result.map_err(|kind| crate::error::record(Error::new(kind, "numerical kernel argument is outside its domain"))).is_ok(); - if let Inputs::Mapped(reads) = inputs && let Some((len, rows)) = by_rows::(reads) { + match inputs.borrowed::(range.clone(), checked) { + Ok(Some(args)) => return output.write(range.len(), &|part, slots, _| unsafe { + slots.fill(part.len(), |out| recorded(apply(args.map(|arg| arg.part(part.start, part.len())), part.len(), out))) + }), + Ok(None) => {}, Err(error) => { output.target.fill(range.len()); return Err(error); } + } + if let Some((len, rows)) = by_rows::(inputs.0) { return output - .write(range.len(), &|part, slots| unsafe { slots.fill(part.len(), |out| recorded(fill_rows(&rows, len, part.start + range.start, out, apply))) }); + .write(range.len(), &|part, slots, _| unsafe { slots.fill(part.len(), |out| recorded(fill_rows(&rows, len, part.start + range.start, out, apply))) }); } - output.write(range.len(), &|part, slots| { + output.write(range.len(), &|part, slots, work| { let part = part.start + range.start..part.end + range.start; - let len = if inputs.in_place::(&part) { part.len().max(1) } else { crate::CHUNK }; - let mut scratch: [Vec; N] = std::array::from_fn(|_| Vec::new()); - part.clone().step_by(len).fold(true, |ok, at| { + let in_place = inputs.in_place::(&part); + let len = if in_place { part.len().max(1) } else { crate::CHUNK }; + let mut local: [Vec; N] = std::array::from_fn(|_| Vec::new()); + let (mut work, mut copied) = (work, 0); + let mut scratch = match work.as_deref_mut() { Some(work) if !in_place => work.inputs::(), _ => local.each_mut() }; + let ok = part.clone().step_by(len).fold(true, |ok, at| { let run = at..(at + len).min(part.end); let n = run.len(); let args = match inputs.part::(run, checked, &mut scratch) { - Ok(args) => args, + Ok((args, n)) => { copied += n; args } Err(error) => { crate::error::record(error); return unsafe { slots.fill(n, |out| { for slot in out { slot.write(O::FILL); } false }) } & ok; } }; (unsafe { slots.fill(n, |out| recorded(apply(args, n, out))) }) & ok - }) + }); + if let Some(work) = work { work.record_copies(copied); } + ok }) } @@ -139,15 +145,21 @@ fn map_blocks( const ROWS: usize = 16; /// The row length `len` and each input as rows of `len` items (`Read::rows`), for a map that reads faster by rows than by blocks: at least two rows of at least `ROWS` items, and an input that blocks would copy, because each of its rows reads one item that changes from row to row, or its rows don't follow on from each other, as an outer product's operands do. `len` is the `k` of the first mapped input's repeat or cycle of `k` items, and else the last axis that input steps along. `None` otherwise. -fn by_rows<'a, F: crate::buffer::Cast>(inputs: &[crate::agreement::Read<'a>]) -> Option<(usize, Vec>)> { - use crate::agreement::{Mapping, ReadRoute}; +fn by_rows<'a, F: crate::buffer::Cast, const N: usize>(inputs: &[crate::agreement::Read<'a>]) -> Option<(usize, [crate::agreement::Rows<'a, F>; N])> { + use crate::{agreement::{Mapping, ReadRoute}, gather::Steps}; let len = inputs.iter().find_map(|read| match read.plan.route() { - ReadRoute::Map { map: Mapping::Linear(k) | Mapping::Tiled(k), .. } => Some(*k), ReadRoute::Map { .. } => read.plan.shape_hint().last().copied(), + ReadRoute::Map(mapped) => match &mapped.map { + Mapping::Linear(k) | Mapping::Tiled(k) => Some(*k), + Mapping::Steps { steps, .. } if Steps::count(steps) == read.plan.len() => steps.last().map(Steps::len), + _ => Some(read.plan.len()), + }, _ => None })?; let items = inputs.first()?.plan.len(); if len < ROWS || items <= len || items % len != 0 { return None; } - let rows = inputs.iter().map(|read| read.rows([items / len, len])).collect::>>()?; + let mut rows = std::array::from_fn(|_| None); + for k in 0..N { rows[k] = Some(inputs.get(k)?.rows([items / len, len])?); } + let rows: [crate::agreement::Rows<'a, F>; N] = rows.map(Option::unwrap); rows.iter().any(|row| if row.run { row.step != len as isize } else { row.step != 0 }).then_some((len, rows)) } @@ -169,9 +181,9 @@ fn fill_rows( } pub(super) fn copied(inputs: Inputs<'_, '_>, range: std::ops::Range, checked: bool, output: &mut Output<'_, '_>) -> Result<(), Error> { - with_dtype!(output.dtype, T => map_values::(inputs, range, checked, false, output, |x| x[0], |_, _| None)) + with_dtype!(output.result_dtype, T => map_values::(inputs, range, checked, false, output, |x| x[0], |_, _| None)) } pub(super) fn constant(value: bool, n: usize, output: &mut Output<'_, '_>) -> Result<(), Error> { - output.write::(n, &|part, slots| slots.pairs(Arg::One(value), Arg::One(()), part.len(), |x, ()| (x, true))) + output.write::(n, &|part, slots, _| slots.pairs(Arg::One(value), Arg::One(()), part.len(), |x, ()| (x, true))) } diff --git a/kernels/src/operation/map/complex.rs b/kernels/src/operation/map/complex.rs index aaa8b41..9a32742 100644 --- a/kernels/src/operation/map/complex.rs +++ b/kernels/src/operation/map/complex.rs @@ -10,7 +10,7 @@ where K: crate::facts::Formula, F: crate::Native + crate::buffer::Cast + Outp map_values::, bool, N>( inputs, range, variant.checked, judge.checked, output, |x| value(x).re != F::zero(), |x, r| invalid(x, num_complex::Complex::new(if r { F::one() } else { F::zero() }, F::zero()))) - } else if output.dtype.is_real() { + } else if output.result_dtype.is_real() { map_values::, F, N>( inputs, range, variant.checked, judge.checked, output, |x| value(x).re, |x, r| invalid(x, num_complex::Complex::new(r, F::zero()))) } else { @@ -21,11 +21,13 @@ where K: crate::facts::Formula, F: crate::Native + crate::buffer::Cast + Outp } } -/// A formula's complex results from real arguments (`Formula::real_complex`). An item that it declines takes `Formula::complex` on its arguments as complex numbers. +/// Complex results from real arguments through a native slice, or `Formula::real_complex` with `Formula::complex` for the items it declines. pub(super) fn map_real_complex( - object: &K, variant: crate::domain::Variant, inputs: Inputs<'_, '_>, range: std::ops::Range, ops: &NativeOps<{ FUSED }>, output: &mut Output<'_, '_> + object: &K, variant: crate::domain::Variant, inputs: Inputs<'_, '_>, range: std::ops::Range, ops: &NativeOps<{ FUSED }>, + slice: Option>>, output: &mut Output<'_, '_> ) -> Result<(), Error> where K: crate::facts::Formula, F: crate::Native + crate::buffer::Cast, num_complex::Complex: OutputItem { map_blocks::, N>(inputs, range, variant.checked, output, &|args, n, out| { + if let Some(slice) = slice { slice(object, ops, args, n, out); return Ok(()); } let nan = num_complex::Complex::new(F::nan(), F::nan()); if !crate::slice::mapped(&args, n, out, #[inline(always)] |x| object.real_complex(ops, x).map_or((nan, false), |z| (z, true))) { for (i, z) in unsafe { crate::slice::written(&mut out[..n]) }.iter_mut().enumerate() { diff --git a/kernels/src/operation/map/io.rs b/kernels/src/operation/map/io.rs index c50d492..1025880 100644 --- a/kernels/src/operation/map/io.rs +++ b/kernels/src/operation/map/io.rs @@ -1,111 +1,34 @@ use super::*; -use crate::{buffer::Cast, parallel::Slots}; +use crate::buffer::Cast; #[derive(Clone, Copy)] -pub(crate) enum Inputs<'a, 'b> { Mapped(&'a [crate::agreement::Read<'b>]), Ready(&'a [View<'b>]) } +pub(crate) struct Inputs<'a, 'b>(pub(crate) &'a [crate::agreement::Read<'b>]); impl<'a, 'b> Inputs<'a, 'b> { - pub(crate) fn len(self) -> usize { match self { Self::Mapped(x) => x.len(), Self::Ready(x) => x.len() } } - pub(crate) fn selected(self, k: usize) -> Self { match self { Self::Mapped(x) => Self::Mapped(&x[k..k + 1]), Self::Ready(x) => Self::Ready(&x[k..k + 1]) } } - pub(crate) fn real(self) -> bool { - let real = |v: View<'_>| !matches!(v.dtype(), DType::Complex(_)); - match self { Self::Mapped(x) => x.iter().all(|r| r.sources.iter().all(|s| real(s.view))), Self::Ready(x) => x.iter().all(|&v| real(v)) } - } - pub(crate) fn in_place(self, range: &std::ops::Range) -> bool { - match self { Self::Mapped(x) => x.iter().all(|r| r.in_place::(range)), Self::Ready(x) => x.iter().all(|&v| F::of(v).is_some()) } - } - pub(crate) fn part<'s, F: Cast, const N: usize>( - self, range: std::ops::Range, checked: bool, scratch: &'s mut [Vec; N]) -> Result<[Arg<'s, F>; N], Error> where 'a: 's, 'b: 's { + pub(crate) fn len(self) -> usize { self.0.len() } + pub(crate) fn selected(self, k: usize) -> Self { Self(&self.0[k..k + 1]) } + pub(crate) fn real(self) -> bool { self.0.iter().all(|r| r.sources.iter().all(|s| !matches!(s.view.dtype(), DType::Complex(_)))) } + pub(crate) fn in_place(self, range: &std::ops::Range) -> bool { self.0.iter().all(|r| r.in_place::(range)) } + pub(crate) fn borrowed(self, range: std::ops::Range, checked: bool) -> Result; N]>, Error> { let mut args = [Arg::Many(&[]); N]; - for (k, (arg, scratch)) in args.iter_mut().zip(scratch).enumerate() { - *arg = match self { - Self::Mapped(x) => x[k].part(range.clone(), checked, scratch)?, - Self::Ready(x) => { - let view = if x[k].len() == 1 { x[k] } else { x[k].slice(range.clone()) }; - match F::of(view) { - Some(items) => Arg::of(items), - None => { - scratch.clear(); - F::extend_cast(view, view.len(), checked, scratch)?; - Arg::of(scratch) - } - } - } - }; - } - Ok(args) - } -} - -/// Where a map appends its results: an owned buffer or the next slots of a caller's result. -pub enum MapOutput<'a, 'b> { - Buffer(&'a mut Buffer), F16(&'a mut Slots<'b, half::f16>), F32(&'a mut Slots<'b, f32>), F64(&'a mut Slots<'b, f64>), Bool(&'a mut Slots<'b, bool>), - U8(&'a mut Slots<'b, u8>), I16(&'a mut Slots<'b, i16>), I32(&'a mut Slots<'b, i32>), I64(&'a mut Slots<'b, i64>), - C32(&'a mut Slots<'b, num_complex::Complex32>), C64(&'a mut Slots<'b, num_complex::Complex64>), Char(&'a mut Slots<'b, char>) -} - -pub(crate) trait OutputItem: Element { fn slots<'s, 'b>(output: &'s mut MapOutput<'_, 'b>) -> Option<&'s mut Slots<'b, Self>>; } - -macro_rules! outputs { - ($($v:ident: $t:ty),+ $(,)?) => { - impl MapOutput<'_, '_> { - fn dtype(&self) -> DType { match self { Self::Buffer(b) => b.dtype(), $(Self::$v(_) => <$t>::DTYPE),+ } } - fn copy(&mut self, view: View<'_>) { - match self { Self::Buffer(b) => { b.extend_from(view, 0..view.len()); }, $(Self::$v(s) => s.copy(<$t>::of(view).expect("result conversion has the destination's dtype"))),+ } - } - pub(crate) fn fill(&mut self, n: usize) { - match self { - Self::Buffer(b) => with_dtype!(b.dtype(), T => { let v = T::vec_mut(b).expect("a buffer has its dtype"); v.resize(v.len() + n, T::FILL); }), - $(Self::$v(s) => { s.pairs(Arg::One(<$t>::FILL), Arg::One(()), n, |x, ()| (x, true)); }),+ - } - } - } - $(impl OutputItem for $t { - fn slots<'s, 'b>(output: &'s mut MapOutput<'_, 'b>) -> Option<&'s mut Slots<'b, Self>> { - match output { MapOutput::$v(s) => Some(s), _ => None } - } + for (arg, read) in args.iter_mut().zip(self.0) { + let Some(items) = read.borrowed(range.clone(), checked)? else { return Ok(None) }; + *arg = items; } - impl<'a, 'b> From<&'a mut Slots<'b, $t>> for MapOutput<'a, 'b> { - fn from(slots: &'a mut Slots<'b, $t>) -> Self { MapOutput::$v(slots) } - })+ - }; -} -outputs!(F16: half::f16, F32: f32, F64: f64, Bool: bool, U8: u8, I16: i16, I32: i32, I64: i64, - C32: num_complex::Complex32, C64: num_complex::Complex64, Char: char); - -pub(crate) struct Output<'a, 'b> { pub target: MapOutput<'a, 'b>, pub dtype: DType, pub checked: bool } - -impl Output<'_, '_> { - /// The block body is already erased here: output conversion never instantiates a formula at a second computing type. - pub(crate) fn write(&mut self, n: usize, apply: &(dyn Fn(std::ops::Range, &mut Slots<'_, O>) -> bool + Sync)) -> Result<(), Error> { - let checked = self.checked; - if let MapOutput::Buffer(buffer) = &mut self.target { - if let Some(values) = O::vec_mut(buffer) { crate::parallel::extended(values, n, crate::CHUNK, apply); } else { - with_dtype!(buffer.dtype(), T => { - let values = T::vec_mut(buffer).expect("a buffer has its dtype"); - crate::parallel::extended(values, n, crate::CHUNK, &|part, slots| { - converted::(&mut slots.into(), part, checked, apply).map_err(crate::error::record).is_ok() - }); - }); - } - return Ok(()); - } - if let Some(slots) = O::slots(&mut self.target) { apply(0..n, slots); return Ok(()); } - converted(&mut self.target, 0..n, checked, apply) + Ok(Some(args)) } -} - -fn converted( - target: &mut MapOutput<'_, '_>, range: std::ops::Range, checked: bool, apply: &(dyn Fn(std::ops::Range, &mut Slots<'_, O>) -> bool + Sync) -) -> Result<(), Error> { - let mut scratch = Vec::::new(); - for at in (range.start..range.end).step_by(crate::CHUNK) { - let n = crate::CHUNK.min(range.end - at); - scratch.clear(); - Slots::append(&mut scratch, n, |slots| apply(at..at + n, slots)); - match Buffer::cast_results(O::view(&scratch), target.dtype(), checked) { - Ok(buffer) => target.copy(buffer.view()), Err(error) => { target.fill(range.end - at); return Err(error); } + pub(crate) fn part<'s, F: Cast, const N: usize>( + self, range: std::ops::Range, checked: bool, scratch: &'s mut [&mut Vec; N] + ) -> Result<([Arg<'s, F>; N], usize), Error> where 'a: 's, 'b: 's { + let (mut args, mut copied) = ([Arg::Many(&[]); N], 0); + for (k, scratch) in scratch.iter_mut().enumerate() { + if let Some(j) = self.0[..k].iter().position(|read| + std::ptr::eq(read.plan, self.0[k].plan) && std::ptr::eq(read.sources, self.0[k].sources) + ) { args[k] = args[j]; continue; } + let (items, n) = self.0[k].part_counted(range.clone(), checked, scratch)?; + args[k] = items; + copied += n; } + Ok((args, copied)) } - Ok(()) } diff --git a/py/src/explain.rs b/py/src/explain.rs index dcf42d2..5f5ce19 100644 --- a/py/src/explain.rs +++ b/py/src/explain.rs @@ -1,6 +1,6 @@ //! Python access to the shared producer-event capture. use pyo3::{prelude::*, types::PyDict}; -use rustymath::events::{self, Event, EventKind, Field, Report, Run, Sink}; +use rustymath::events::{Event, EventKind, Field, Report, Run, Sink}; use std::sync::{Arc, Mutex, MutexGuard, PoisonError}; #[derive(Default)] @@ -30,9 +30,9 @@ impl Explanation { report.records.iter().map(|event| { let record = PyDict::new(py); record.set_item("kind", event.name())?; - for ((name, _), field) in events::TABLES[event.table()].1.iter().zip(event.fields()) { + for ((name, _), field) in event.table().1.iter().zip(event.fields()) { match field { - Field::Text(t) => record.set_item(name, t)?, Field::Count(n) => record.set_item(name, n)?, Field::Truth(t) => record.set_item(name, t)? + Field::Text(t) => record.set_item(name, t.as_ref())?, Field::Count(n) => record.set_item(name, n)?, Field::Truth(t) => record.set_item(name, t)? } } Ok(record.unbind()) diff --git a/pyproject.toml b/pyproject.toml index 88ab320..52dfa76 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -13,7 +13,7 @@ authors = [{name = "Jeremy Howard", email = "github@jhoward.fastmail.fm"}] classifiers = ["Programming Language :: Rust", "Programming Language :: Python :: Implementation :: CPython"] [project.optional-dependencies] -dev = ["fastship>=0.0.11", "fastws-cli>=0.0.34", "maturin>=1.0,<2.0", "pytest", "pytest-timeout", "nbdev", "numpy", "pandas", "scipy"] +dev = ["fastship>=0.0.11", "fastws-cli>=0.0.35", "maturin>=1.0,<2.0", "pytest", "pytest-timeout", "nbdev", "numpy", "pandas", "scipy"] [project.urls] Homepage = "https://github.com/AnswerDotAI/rustymath" @@ -36,7 +36,7 @@ module-name = "rustymath._core" macos-deployment-target = "15.0" [tool.uv] -cache-keys = [{ file = "pyproject.toml" }, { file = "Cargo.toml" }, { file = "core/Cargo.toml" }, { file = "kernels/Cargo.toml" }, { file = "py/Cargo.toml" }] +cache-keys = [{ file = "pyproject.toml" }, { file = "Cargo.toml" }, { file = "core/Cargo.toml" }, { file = "kernels/Cargo.toml" }, { file = "py/Cargo.toml" }, { file = "timing/Cargo.toml" }, { file = "buffer/Cargo.toml" }] [tool.pytest.ini_options] testpaths = ["tests"] diff --git a/src/create.rs b/src/create.rs index c9276e3..93611ec 100644 --- a/src/create.rs +++ b/src/create.rs @@ -85,6 +85,6 @@ pub fn linspace(start: f64, stop: f64, n: usize, endpoint: bool) -> Result Result, Error> { let exponents = linspace(start, stop, n, endpoint)?; let mut out = Vec::with_capacity(n); - crate::pairs(crate::Arg::One(base), crate::Arg::Many(&exponents), n, &mut out, |x, y| (crate::facts::formula_at(&crate::Pow, [x, y]), true)); + rustymath_core::slice::pairs(crate::Arg::One(base), crate::Arg::Many(&exponents), n, &mut out, |x, y| (crate::facts::formula_at(&crate::Pow, [x, y]), true)); Ok(out) } diff --git a/src/lib.rs b/src/lib.rs index 53d4222..2c16d2d 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -2,7 +2,31 @@ #![doc = include_str!("lib.md")] #[macro_use] extern crate rustymath_core; -pub use rustymath_core::*; +// The public surface is math, arrays and their execution facts, not every helper +// made public inside core so the kernels crate can call it. +pub use rustymath_core::{ + agreement, buffer, characters, convert, derivative, distribution, domain, error, events, exact, facts, gather, hashing, integer, interrupt, + layout, memory, number_theory, optimize, parallel, polynomial, random, recipe, reduce, search, signal, stats, structure, tolerant, + fused, unfused, half, num_complex, +}; +pub use rustymath_core::{ + FUSED, Float, Native, Arg, Unary, Binary, Ternary, + Abs, Acos, Acosh, Add, Angle, Asin, Asinh, Atan, Atan2, Atanh, Betainc, BetaincAt, Betaincinv, Binom, BinomDk, BinomDn, + BitAnd, BitOr, Bitwise, Cbrt, Choose, Cis, Clip, Conj, Cos, Cosh, Digamma, Divide, Erf, Erfc, Erfcinv, Erfinv, Exp, Expm1, + Factorial, Gamma, Gammainc, Gammaincc, Gammaincinv, Hypot, Imag, Isfinite, Isinf, Isnan, LambertW, LeftShift, Lgamma, Ln, Log1p, Log2, + Log10, Logb, MathishDivide, MathishLogb, Multiply, Nand, Negate, Nor, Not, PiTimes, Polygamma, Pow, Real, RealAcos, RealAcosh, RealAsin, + RealAtanh, RealLn, RealPow, RealSqrt, RealSqrt1mx2, RealSqrtm1mx2, RealSqrtx2m1, Reciprocal, RightShift, Rint, Root, Round, Scale, + Sign, Sin, Sinh, Sqrt, Sqrt1mx2, Sqrt1px2, Sqrtm1mx2, Sqrtx2m1, Subtract, Tan, Tanh, Trunc, Exponent, +}; +pub use rustymath_core::{ + acos, acosh, asin, asinh, atan, atan2, atanh, beta, betainc, betaincinv, binom, binom_dk, binom_dn, bitwise, cbrt, cos, cosh, digamma, + divide, erf, erfc, erfcinv, erfinv, exp, expm1, factorial, gamma, gammainc, gammaincc, gammaincinv, hypot, lambertw, lgamma, ln, log1p, + log2, log10, logb, polygamma, pow, root, sin, sinh, sqrt1mx2, sqrt1px2, sqrtx2m1, tan, tanh, +}; +pub use rustymath_core::{ + element_types, function_doc, native_reduction, with_complex_width, with_complexes, with_dtype, with_float_width, with_floats, with_inexact, + with_int_pair, with_ints, with_reduce, with_tolerance, with_width, +}; pub use rustymath_kernels::operation; pub mod array; pub mod create; diff --git a/tests/accuracy.rs b/tests/accuracy.rs index 8a7210a..c1ad28e 100644 --- a/tests/accuracy.rs +++ b/tests/accuracy.rs @@ -686,7 +686,7 @@ mod complex { cases.sort_by(|a, b| a.0.cmp(&b.0)); for group in cases.chunk_by(|a, b| a.0 == b.0) { let name = &group[0].0; - let count = if bulk && name == "exp" { group.len().max(rustymath::CHUNK + 1) } else { group.len() }; + let count = if bulk && name == "exp" { group.len().max(rustymath_core::CHUNK + 1) } else { group.len() }; let group = group.iter().cycle().take(count).collect::>(); let inputs = (0..group[0].1.len()).map(|i| group.iter().map(|(_, a, _)| c(narrow(a[i].re), narrow(a[i].im))).collect()).collect::>(); let got = apply_many(name, &inputs, bulk).expect("a reference function"); @@ -702,10 +702,10 @@ mod complex { out } - fn bulk_edges() where Complex: Element { + fn bulk_edges() where Complex: Element { let real = |x| c(R::narrow(x), R::zero()); assert!(evaluate::(Exp, [&[]], true).is_empty()); - let n = rustymath::CHUNK + 1; + let n = rustymath_core::CHUNK + 1; let bases = [real(0.0), real(2.0), real(-2.0)].into_iter().cycle().take(n).collect::>(); for (exponent, values) in [(0.0, [1.0, 1.0, 1.0]), (1.0, [0.0, 2.0, -2.0]), (2.0, [0.0, 4.0, 4.0]), (-1.0, [f64::INFINITY, 0.5, -0.5])] { let got = evaluate(Pow, [&bases, &[real(exponent)]], true); @@ -714,6 +714,18 @@ mod complex { let exponents = [0.0, 1.0, 2.0, -1.0].into_iter().cycle().take(n).map(real).collect::>(); let got = evaluate(Pow, [&[real(-2.0)], &exponents], true); assert!(got.iter().zip([1.0, -2.0, 4.0, -0.5].into_iter().cycle()).all(|(&g, w)| g == real(w)), "broadcast base"); + let cis = Kernel::resolve(rustymath::Cis, Variant::checked(Complex::::DTYPE), [R::DTYPE]).unwrap(); + let ys = [0.0, FRAC_PI_2, -PI, 0.3, 1e20].map(R::narrow).into_iter().cycle().take(n).collect::>(); + let tolerance = 8.0 * if R::WIDE { f64::EPSILON } else { f32::EPSILON.into() }; + for (count, xs) in [(0, &ys[..0]), (n, &ys[4..5]), (n, &ys[..])] { + let got: Vec> = crate::support::mapped(&cis, [xs], count).unwrap(); + assert_eq!(got.len(), count); + for (i, z) in got.iter().enumerate() { + let want = unary(rustymath::Cis, c(xs[if xs.len() == 1 { 0 } else { i }], R::zero())); + assert!((z.re.wide() - want.re.wide()).hypot(z.im.wide() - want.im.wide()) <= tolerance); + } + if xs.len() > 1 { assert!(got[..3] == [real(1.0), c(R::zero(), R::one()), real(-1.0)]); } + } } #[test] diff --git a/tests/linalg.rs b/tests/linalg.rs index aca0586..2de7919 100644 --- a/tests/linalg.rs +++ b/tests/linalg.rs @@ -39,7 +39,7 @@ mod products { let fk = Kernel::resolve(f, Variant::checked(gk.result_dtype()), [gk.result_dtype(); 2]).unwrap(); inner_with( a, b, $size, |x, y, agreement| { - assert!(agreement.len <= rustymath::CHUNK); + assert!(agreement.len <= rustymath_core::CHUNK); let sources = [Binding::full(x), Binding::full(y)]; let plans = [ ReadPlan::input(0, x.len()).remapped(agreement.left.clone(), agreement.len, Missing::Default), @@ -53,7 +53,7 @@ mod products { ) }}; } - let [m, k, n] = [2, rustymath::CHUNK + 3, 3]; + let [m, k, n] = [2, rustymath_core::CHUNK + 3, 3]; let a = (0..m * k).map(|i| (i % 7) as i64 - 3).collect::>(); let b = (0..k * n).map(|i| (i % 5) as i64 - 2).collect::>(); let out = inner!(rustymath::Subtract, rustymath::Multiply, i64::view(&a), i64::view(&b), [m, k, n]).unwrap(); @@ -68,8 +68,8 @@ mod products { (0..k - 2).rev().fold(t, |t, p| term(p) == i64::from(t)) }).collect::>(); assert_eq!(bool::of(out.view()).unwrap(), expected); - let wide = vec![1i64; 3 * (rustymath::CHUNK + 1)]; - let out = inner!(rustymath::Subtract, rustymath::Multiply, i64::view(&[1, 2, 3]), i64::view(&wide), [1, 3, rustymath::CHUNK + 1]).unwrap(); + let wide = vec![1i64; 3 * (rustymath_core::CHUNK + 1)]; + let out = inner!(rustymath::Subtract, rustymath::Multiply, i64::view(&[1, 2, 3]), i64::view(&wide), [1, 3, rustymath_core::CHUNK + 1]).unwrap(); assert!(i64::of(out.view()).unwrap().iter().all(|&x| x == 2)); let out = inner!(rustymath::Add, Greater(Tolerance::EXACT), i64::view(&[2, 3, 4]), i64::view(&[1, 5, 2]), [1, 3, 1]).unwrap(); assert_eq!(i64::of(out.view()).unwrap(), [2]); @@ -241,7 +241,7 @@ mod products { #[test] fn wrapped_and_half_products_match_folds() { // `Wrapped(PlusTimes)` folds with wrapping `i32` steps whatever its items' type, whether they compute as exact floats or on the `i32` kernel: `u8`, `i16` and `i32` items give the folds of their `i32` values, and the large `i16` and `i32` items overflow 32 bits. Small and moderate `i32` items take the exact float route, and the moderate ones still wrap at 32 bits. f16 `+.×` gives the f32 folds of the same values, which small whole numbers keep exact. - use rustymath::Wrapped; + use rustymath_core::Wrapped; fn widened>(x: &[T]) -> Vec { x.iter().map(|&v| v.into()).collect() } for size @ [m, k, n] in [[0, 3, 4], [3, 0, 4], [3, 4, 0], [5, 7, 9], [MC + 3, KC + 5, 37]] { let wrapped = |a: &[i32], b: &[i32]| folded(a, b, size, 0i32, |s, x, y| s.wrapping_add(x.wrapping_mul(y))); diff --git a/tests/numeric.rs b/tests/numeric.rs index cca2549..4c94417 100644 --- a/tests/numeric.rs +++ b/tests/numeric.rs @@ -575,29 +575,49 @@ mod operation { } #[test] - fn ready_maps_append_direct_converted_and_failed_results() { - use rustymath::{Add, Sqrt, operation::MapOutput, parallel::Slots, tolerant::{Less, Tolerance}}; - let settings = Settings { tolerance: Tolerance::EXACT }; + fn mapped_kernels_append_direct_converted_and_failed_results() { + use rustymath::{Add, Sqrt, buffer::{MapOutput, Scratch}, parallel::Slots, tolerant::{Less, Tolerance}}; + fn append(kernel: &Kernel, inputs: &[rustymath::buffer::View<'_>], n: usize, values: &mut Vec) -> Result<(), Error> { + let layout = Layout::from([n]); + let inputs: Vec<_> = inputs.iter().map(|&data| view(&layout, data, &if data.len() == 1 { Mapping::Single } else { Mapping::Linear(1) })).collect(); + let reads: Vec<_> = inputs.iter().map(ArrayView::read).collect(); + let map = rustymath::operation::Map::Kernel(kernel.clone(), Settings { tolerance: Tolerance::EXACT }); + Slots::append(values, n, |slots| map.run(&reads, 0..n, &mut Scratch::default(), MapOutput::F64(slots))) + } let float = f64::DTYPE; let add = Kernel::resolve(Add, Variant::checked(float), [float; 2]).unwrap(); + let row_shape = Layout::from([2, 32]); + let rows = view(&row_shape, f64::view(&[1., 2.]), &Mapping::Linear(32)); + assert!(add.map(&[rows.read(), rows.read()], 0..65, Settings { tolerance: Tolerance::EXACT }, &mut Buffer::with_capacity(float, 0)).is_err()); let mut values = vec![10.]; - Slots::append(&mut values, 3, |slots| add.map_ready(&[f64::view(&[1.]), f64::view(&[2., 3., 4.])], 3, settings, MapOutput::F64(slots))).unwrap(); + append(&add, &[f64::view(&[1.]), f64::view(&[2., 3., 4.])], 3, &mut values).unwrap(); assert_eq!(values, [10., 3., 4., 5.]); let less = Kernel::resolve(Less(Tolerance::EXACT), Variant::checked(float), [float; 2]).unwrap(); - Slots::append(&mut values, 3, |slots| less.map_ready(&[f64::view(&[3.]), f64::view(&[2., 3., 4.])], 3, settings, MapOutput::F64(slots))).unwrap(); + append(&less, &[f64::view(&[3.]), f64::view(&[2., 3., 4.])], 3, &mut values).unwrap(); assert_eq!(values, [10., 3., 4., 5., 0., 0., 1.]); let complex = num_complex::Complex64::DTYPE; let sqrt = Kernel::resolve(Sqrt, Variant::checked(complex), [complex]).unwrap(); - let result = Slots::append(&mut values, 1, |slots| { - sqrt.map_ready(&[num_complex::Complex64::view(&[num_complex::Complex64::new(-1., 0.)])], 1, settings, MapOutput::F64(slots)) - }); + let result = append(&sqrt, &[num_complex::Complex64::view(&[num_complex::Complex64::new(-1., 0.)])], 1, &mut values); assert!(result.is_err()); assert_eq!(values.len(), 8); let ints = Kernel::resolve(Add, Variant::checked(i64::DTYPE), [float; 2]).unwrap(); let inputs = [f64::view(&[f64::NAN]), f64::view(&[1.])]; - Slots::append(&mut values, 0, |slots| ints.map_ready(&inputs, 0, settings, MapOutput::F64(slots))).unwrap(); - assert!(Slots::append(&mut values, 1, |slots| ints.map_ready(&inputs, 1, settings, MapOutput::F64(slots))).is_err()); + append(&ints, &inputs, 0, &mut values).unwrap(); + assert!(append(&ints, &inputs, 1, &mut values).is_err()); assert_eq!(values.len(), 9); + let data = f64::view(&[-1.5, 2.9]); + assert_eq!(i16::vec(Buffer::floored(data, i16::DTYPE, Tolerance::EXACT).unwrap()).unwrap(), [-2, 2]); + assert_eq!(i16::vec(Buffer::cast_results(data, i16::DTYPE, true).unwrap()).unwrap(), [-1, 2]); + assert_eq!(u8::vec(Buffer::floored(i64::view(&[300]), u8::DTYPE, Tolerance::EXACT).unwrap()).unwrap(), [44]); + let errors = rustymath::error::ErrorSlot::default(); + errors.within(|| { + assert!(Buffer::kept(data, i16::DTYPE).is_ok()); + assert_eq!(errors.take().unwrap().kind, rustymath::error::Kind::Domain); + assert!(Buffer::cast_results(i64::view(&[300]), u8::DTYPE, true).is_ok()); + assert_eq!(errors.take().unwrap().kind, rustymath::error::Kind::Overflow); + assert!(Buffer::cast_results(i64::view(&[0xd800]), char::DTYPE, true).is_ok()); + assert_eq!(errors.take().unwrap().kind, rustymath::error::Kind::Domain); + }); } #[test] @@ -613,6 +633,8 @@ mod operation { let mut out = Buffer::with_capacity(add.result_dtype(), 3); crate::support::map(&add, [bool::view(&[true]), bool::view(&[true, false, true])], 3, &mut out).unwrap(); assert_eq!(i64::vec(out).unwrap(), [2, 1, 2]); + let byte = Kernel::resolve(Add, Variant::checked(u8::DTYPE), [u8::DTYPE; 2]).unwrap(); + assert!(byte.range_on(&[Range::of(Some(200.)); 2], false).hi >= 400.); let checked = Kernel::resolve(Add, Variant::checked(wide), [wide; 2]).unwrap(); let slot = ErrorSlot::new(); let mut out = Buffer::with_capacity(wide, 1); @@ -754,7 +776,7 @@ mod operation { assert_eq!(i64::vec(Buffer::cast_results(out.view(), dtype, true).unwrap()).unwrap(), [expected]); } for inner in [1, 2] { - let outer = rustymath::CHUNK / 3 + 1; + let outer = rustymath_core::CHUNK / 3 + 1; let data: Vec = (0..outer * 3 * inner).map(|i| [3, 2, 1][i / inner % 3]).collect(); let input = [Binding::full(i64::view(&data))]; let plan = ReadPlan::input(0, data.len()); @@ -816,12 +838,31 @@ mod operation { #[test] fn logical_views_compose_prepare_and_run() { let (shape, data) = (Layout::from([3, 2]), [0., 1., 2., 3., 4., 5.]); + let extended = Layout::from([8]); + let extended = view(&extended, f64::view(&data), &Mapping::Linear(1)); + assert_eq!(extended.read().plan.contiguous(1..4), Some((0, 1..4))); + assert_eq!(extended.read().plan.contiguous(1..1), None); + let fallback = ReadPlan::copy(1, 0, Mapping::Linear(1), 2, Missing::Input(ReadPlan::input(0, 1))); + assert!(fallback.contiguous(0..2).is_none()); + assert_eq!(extended.items::().unwrap().as_ref(), &[0., 1., 2., 3., 4., 5., 0., 0.]); + assert!(matches!(extended.read().part_counted::(2..3, true, &mut Vec::new()), Ok((rustymath::Arg::One(2), 0)))); + let sources = [Binding::window(f64::view(&data[2..4]), 2, data.len()).unwrap()]; + let source = ReadPlan::input(0, data.len()); + assert!(Read { sources: &sources, plan: &source }.rows::([3, 2]).is_none()); + let slice = source.remapped(Mapping::Steps { base: 2, steps: vec![Steps::along(1, 0), Steps::along(2, 1)].into() }, 2, Missing::Invalid); + assert_eq!(Read { sources: &sources, plan: &slice }.rows::([1, 2]).unwrap().items, &data[2..4]); + assert!(Read { sources: &sources, plan: &slice }.rows::([1, 1]).is_none()); + let broadcast = Mapping::Steps { base: 0, steps: vec![Steps::along(2, 2), Steps::along(3, 0), Steps::along(2, 1)].into() }; + assert_eq!(broadcast.affine(&[2, 3, 2]), Some((0, vec![2, 0, 1]))); + assert!(broadcast.affine(&[6, 2]).is_none()); let transpose = Mapping::Steps { base: 0, steps: vec![Steps::along(3, 1), Steps::along(2, 3)].into() }; let x = view(&shape, f64::view(&data), &transpose); assert!(x.contiguous().is_none()); assert_eq!(x.items::().unwrap().as_ref(), &[0., 3., 1., 4., 2., 5.]); - assert_eq!(transpose.affine(shape.shape()), Some((0, vec![1, 3]))); + let reshaped = Layout::from([1, 3, 2]); + assert_eq!(x.with_layout(&reshaped).unwrap().affine().map(|(_, base, strides)| (base, strides)), Some((0, vec![0, 1, 3]))); let reversed = Mapping::Steps { base: 0, steps: vec![Steps::Stride { start: 4, len: 3, step: -2 }, Steps::along(2, 1)].into() }; + assert_eq!(reversed.affine(shape.shape()), Some((4, vec![-2, 1]))); let composed = transpose.compose(&reversed, shape.shape(), shape.shape()).unwrap(); assert_eq!(view(&shape, f64::view(&data), &composed).items::().unwrap().as_ref(), &[2., 5., 1., 4., 0., 3.]); let add = Kernel::resolve(rustymath::Add, Variant::checked(f64::DTYPE), [f64::DTYPE; 2]).unwrap(); @@ -848,6 +889,7 @@ mod operation { let prefix = Layout::from([2]); let prefix = view(&prefix, f64::view(&data), &strided); assert_eq!(prefix.items::().unwrap().as_ref(), &[3., 2.]); + assert_eq!(prefix.read().plan.contiguous(0..1), Some((0, 4..5))); assert_eq!(prefix.items::().unwrap().as_ref(), &[3, 2]); assert!(view(&Layout::from([0, 3]), f64::view(&data), &transpose).items::().unwrap().is_empty()); assert!(Mapping::Linear(2).affine(&[3]).is_none()); @@ -867,6 +909,12 @@ mod operation { let chosen = chosen.selected(padding, Some(std::sync::Arc::new(f64::buffer(vec![8.])))).unwrap(); let args = [chosen.view(), ArrayView::compact(&Layout::EMPTY, f64::view(&[1.])), ArrayView::compact(&Layout::EMPTY, f64::view(&[10.]))]; assert_eq!(f64::vec(apply(&clip, args).unwrap()).unwrap(), [8., 4., 2., 4., 8.]); + for (p, run) in [(0, 0), (0, 8), (1, 0), (1, 8), (3, 2), (3, 3), (3, 8)] { + let mut got = Vec::new(); + let tiled = rustymath::agreement::tiled_cycle(&[0u8, 1, 2][..p], &Mapping::Tiled(p), run, &mut got); + assert!(got.iter().enumerate().all(|(i, &v)| v == (i % p) as u8)); + assert_eq!(rustymath::agreement::Period::Cycles(p).tiling(run).map(|(_, n)| n), tiled.map(|_| got.len())); + } } #[test] diff --git a/tests/slices.rs b/tests/slices.rs index a8a223d..15c443c 100644 --- a/tests/slices.rs +++ b/tests/slices.rs @@ -131,11 +131,12 @@ mod forms { Lgamma, Erf, Erfc, Erfinv, Erfcinv, LambertW, Sqrt, Sqrt1mx2, Sqrt1px2, Sqrtx2m1, Abs, Sign, Angle; binary: Pow, Root, Hypot, Atan2, Gammainc, Gammaincc, Gammaincinv, Binom, Logb, MathishLogb, Maximum, Minimum; ternary: Betainc, Betaincinv ); - // `Erfcinv`'s f64 slice form places each block's items by region. Several blocks of uniform items, in random order, reach every region. + // The f64 inverse-error slice forms place each block's items by region. Several blocks of uniform items, in random order, reach every region. use rand::RngExt; let mut rng = crate::support::seeded(1); let q: Vec = (0..2000).map(|_| 2.0 * rng.random::()).collect(); unary(Erfcinv, "Erfcinv", &q); + unary(Erfinv, "Erfinv", &q.iter().map(|&v| v - 1.0).collect::>()); // The logarithms' slice forms take each block of 64 items whose every item is positive, normal and finite through a form with no special values. Blocks of such items, some holding a special value, reach both forms. let mut x: Vec = (0..2000).map(|_| 1000.0 * rng.random::()).collect(); for (i, v) in [0.0, -1.0, 1e-310, f64::INFINITY, f64::NAN, 1.0].into_iter().enumerate() { x[300 * i + 5] = v } @@ -867,6 +868,12 @@ mod tolerant { if rustymath::FUSED { assert_ne!(Exactly.residue(0.1, 1.0), 0.0); } assert_eq!(Exactly.residue(1e300, 1e-300), 1e-300); assert!((T.gcd(0.6, 0.4) - 0.2).abs() < 1e-15 && Exactly.gcd(0.6, 0.4) < 1e-15); + let tiny = f64::from_bits(1); + for (x, y, g) in [ + (0.75, 0.5, 0.25), (1.0, 0.75, 0.25), (3.0 * tiny, 9.0 * tiny, 3.0 * tiny), + (f64::MAX, tiny, tiny), (f64::MAX, f64::MAX, f64::MAX), + ] { assert_eq!(Exactly.gcd(x, y), g); } + assert_eq!(Exactly.gcd(0.6, 0.4), f64::EPSILON / 2.0); assert_eq!((T.gcd(12.0, 18.0), T.gcd(-12.0, 18.0), T.gcd(0.0, 2.5), T.gcd(0.0, 0.0), T.gcd(12.0f32, 18.0)), (6.0, 6.0, 2.5, 0.0, 6.0)); assert!(T.gcd(f64::NAN, 1.0).is_nan() && T.gcd(f64::NAN, f64::NAN).is_nan() && T.gcd(f64::INFINITY, 1.0).is_nan()); } @@ -1162,6 +1169,19 @@ mod gather { assert!(check(&[20], &[2, 19], vec![along(2, 1), along(19, 1)])); assert!(check(&[5, 6, 2], &[3, 4, 3, 3, 2], vec![along(3, 12), along(4, 2), along(3, 12), along(3, 2), along(2, 1)])); } + /// Shared index storage still distinguishes its axis interpretation, scale and selected range. + #[test] + fn index_tables_compare_source_and_offsets() { + use rustymath::{buffer::{Buffer, IntBuf}, gather::{IndexSource, Table}}; + let values: std::sync::Arc = std::sync::Arc::new(Buffer::Integers(IntBuf::I16(vec![-1, 0, 1]))); + let table = Table::indices(values.clone(), 4, None); + assert_eq!(table, table.clone()); + assert_eq!(table, Table::new(vec![3, 0, 1])); + assert_ne!(table, table.clone().scaled(2)); + assert_ne!(table, Table::indices(values, 5, None)); + assert_ne!(table.slice(0, 2), table.slice(1, 2)); + } + /// Index arrays pair by broadcasting and take the place of the first one, while a range keeps its axis in place. With one index array, pairing is the cross product. #[test] fn index_arrays_pair() { diff --git a/tests/speed_python.txt b/tests/speed_python.txt index 6d702c4..52c13ea 100644 --- a/tests/speed_python.txt +++ b/tests/speed_python.txt @@ -12,6 +12,13 @@ x*2 2.2 2.5 instructions linux 1.019 -x 2.2 2 instructions linux 0.941 # 7.1 copies, for Exp's 7 copies, plus a tenth of a copy for the Python call. rm.exp(x) 15 9.7 instructions linux 10.44 +# One real input and one complex output per item; bulk cost is uncalibrated. +rm.cis(x) - - +# Error-function slice forms at f64; bulk costs are uncalibrated. +rm.erfc(x) - - +rm.erfinv(x) - - +# Sorted inputs exercise whole-region Erfc blocks at both widths; costs are uncalibrated. +rm.erfc(xs) - - # 1.9 copies, for reading the Boolean mask, x and y, a select, one write, plus a tenth of a copy for the Python call. rm.where(mask, x, y) 3.8 5.6 instructions linux 2.173 # 1.6 copies, for one read, two comparisons, one write, plus a tenth of a copy for the Python call. @@ -72,6 +79,12 @@ np.asarray(x) 0.5 0.5 instructions linux 0.29 xf+yf 3.4 3.2 instructions linux 0.735 # 5.7 copies, for Exp's 5.5 copies at f32, plus the call. rm.exp(xf) 12 6.5 instructions linux 4.193 +# Real-input cis with complex output at f32; bulk cost is uncalibrated. +rm.cis(xf) - - +# Error-function slice forms at f32; bulk costs are uncalibrated. +rm.erfc(xf) - - +rm.erfinv(xf) - - +rm.erfc(xsf) - - # 2.1 copies, for reading the mask, xf and yf, a select, one write,, plus the call. rm.where(mask, xf, yf) 4.2 5.5 instructions linux 1.116 # 0.7 copies, for one read per item, plus the call. diff --git a/tests/test_speed.py b/tests/test_speed.py index 37725c5..231c8ca 100644 --- a/tests/test_speed.py +++ b/tests/test_speed.py @@ -12,7 +12,9 @@ def arrays(): rng = np.random.default_rng(0) u = lambda *shape: rng.uniform(-1, 1, shape) x, y = u(N), u(N) + xs = np.sort(x) return dict(np=np, rm=rm, xn=x, x=rm.array(x), y=rm.array(y), xf=rm.array(x.astype(np.float32)), yf=rm.array(y.astype(np.float32)), + xs=rm.array(xs), xsf=rm.array(xs.astype(np.float32)), m=rm.array(u(256, 256)), a=rm.array(u(256, 256)), b=rm.array(u(256, 256)), row=rm.array(u(256)), idx=rm.array(rng.permutation(N)), img=rm.array(u(256, 256)), k=rm.array(u(3, 3)), kern=rm.array(u(16)), frames=rm.array(u(1024, 64)), r=rm.array(rng.integers(0, 256, N)), c=rm.array(rng.integers(0, 256, N)), mask=rm.array(u(N) > 0)) @@ -46,16 +48,22 @@ def attempt(case, namespace, units, timing): return timing.line(section, code, copies, '-' if budget is None else str(budget), best / N, footprint, (per_item, None)), copies, per_item def retimed(over, namespace, units, timing): - "Time each of `over`'s `(case, line)` pairs again, printing its new line, and give those still over budget, with their new lines." - again = [] - for case, _ in over: - line, copies, _ = attempt(case, namespace, units, timing) - print(line, file=sys.stderr) - if case[2].over(copies): again.append((case, line)) - return again + "Retry cases over budget after each of `timing.pauses`, warming before each retry, and return those still over. Print each failure before its retry and each new measurement." + for pause in timing.pauses: + if not over: break + for _, line in over: print(f'{line}: over, timed again after {pause:g} s', file=sys.stderr) + time.sleep(pause) + timing.warm() + again = [] + for case, _ in over: + line, copies, _ = attempt(case, namespace, units, timing) + print(line, file=sys.stderr) + if case[2].over(copies): again.append((case, line)) + over = again + return over def test_speed(capsys): - "Every case costs no more than its budget, and retires no more instructions than its recorded count allows, on the CPU that the budgets were set for (`timing.checks_budgets`). On any other CPU, it prints every cost, then each case still over, marked `over`, and fails on none. An ordinary run holds the timing lock shared, on one thread, beside builds and other tests. It then times each case over its budget again if it can hold the lock exclusively at once, and otherwise marks the case `unconfirmed` and doesn't fail it. A calibration run, which `RUSTYMATH_COSTS` asks for, holds the lock exclusively, prints every cost and times each case over its budget again after each of `timing.pauses`. Each case whose gap's hold has outlived the gap is printed, marked `stale`, and each whose instructions fall below its record divided by the margin is marked `fewer instructions`." + "Every case costs no more than its budget, and retires no more instructions than its recorded count allows, on the CPU that the budgets were set for (`timing.checks_budgets`). On any other CPU, it prints every cost, then each case still over, marked `over`, and fails on none. An ordinary run holds the timing lock shared, on one thread, beside builds and other tests. It retries each case over its budget after each of `timing.pauses` if it can hold the lock exclusively at once, and otherwise marks the case `unconfirmed` and doesn't fail it. A calibration run, which `RUSTYMATH_COSTS` asks for, holds the lock exclusively, prints every cost and times each case over its budget again after each of `timing.pauses`. Each case whose gap's hold has outlived the gap is printed, marked `stale`, and each whose instructions fall below its record divided by the margin is marked `fewer instructions`." timing = rm.timing(shared=not CALIBRATING) namespace, cases = arrays(), timing.budgets(Path(__file__).with_name('speed_python.txt').read_text()) show, failures = not timing.checks_budgets or CALIBRATING, [] @@ -73,19 +81,13 @@ def test_speed(capsys): if record.over(instructions, timing.count_margin): failures.append(f'{line}: instructions over the {record.value} recorded') elif instructions * timing.count_margin < record.value: print(f'{line}: fewer instructions than the {record.value} recorded', file=sys.stderr) over = [(case, line) for case, line, copies, _ in measured if case[2] is not None and case[2].over(copies)] - for pause in timing.pauses if CALIBRATING else (): - if not over: break - for _, line in over: print(f'{line}: over, timed again after {pause:g} s', file=sys.stderr) - time.sleep(pause) - timing.warm() - over = retimed(over, namespace, units, timing) + if CALIBRATING: over = retimed(over, namespace, units, timing) if over and not CALIBRATING: if (quiet := rm.try_timing()) is None: for _, line in over: print(f'{line}: unconfirmed, because builds or tests hold the timing lock', file=sys.stderr) over = [] else: with held(quiet): - timing.warm() for _, line in over: print(f'{line}: over, timed again with the timing lock held exclusively', file=sys.stderr) over = retimed(over, namespace, units, timing) failures += [line for _, line in over]