sciencekit

Why does sciencekit exist?

Python ML libraries carry decades of ecosystem — along with the GIL, insecure pickle, and layers of interpretation between your code and the metal. sciencekit brings the same algorithmic reach as scikit-learn to a compiled, native base — no garbage collector, no ergonomic compromises:

The API you will write

Mandatory builder pattern, zero-copy inputs, and full descriptive names — no max_iter anywhere:

// Illustrative — requires the sciencekit crates (Phase 0+).
let model = SKKMeansClassifierBuilder::new()
    .number_of_clusters(8)
    .maximum_iterations(300)
    .execution_mode(SKExecutionMode::Automatic)
    .build()?;

model.fit(&training_data_view)?;
let predictions = model.predict(&test_data_view)?;

Inputs take ArrayView/CowArray (or sparse views via sprs) — never copies of huge matrices. Pipelines are validated at compile time.

Try it in your browser

The example below mirrors the exact shape of the sciencekit API — mandatory builder, full descriptive names, execution_mode on every builder. Edit the code and hit Run: it executes for real, right on the page.

// The sciencekit API follows the builder pattern everywhere.
// This standalone demo mirrors its shape — the real crates arrive with Phase 0.
// Bounded options are enums, never strings, so `execution_mode` cannot be
// misspelled.
#[derive(Debug, Clone, Copy)]
enum SKExecutionMode {
    Automatic,
    InProcessSynchronous,
    OutOfCoreStreaming,
}

#[derive(Debug)]
struct TrainedKMeansModel {
    number_of_clusters: usize,
    maximum_iterations: usize,
    execution_mode: SKExecutionMode,
}

struct KMeansBuilder {
    number_of_clusters: usize,
    maximum_iterations: usize,
    execution_mode: SKExecutionMode,
}

impl KMeansBuilder {
    fn new() -> Self {
        Self {
            number_of_clusters: 3,
            maximum_iterations: 100,
            execution_mode: SKExecutionMode::Automatic,
        }
    }

    fn number_of_clusters(mut self, value: usize) -> Self {
        self.number_of_clusters = value;
        self
    }

    fn maximum_iterations(mut self, value: usize) -> Self {
        self.maximum_iterations = value;
        self
    }

    fn execution_mode(self, mode: SKExecutionMode) -> Self {
        println!("execution mode requested: {mode:?}");
        self
    }

    fn build(self) -> TrainedKMeansModel {
        TrainedKMeansModel {
            number_of_clusters: self.number_of_clusters,
            maximum_iterations: self.maximum_iterations,
            execution_mode: self.execution_mode,
        }
    }
}

fn main() {
    let model = KMeansBuilder::new()
        .number_of_clusters(8)
        .maximum_iterations(300)
        .execution_mode(SKExecutionMode::Automatic)
        .build();

    println!("model ready: {model:?}");
}

Who is it for?

• ML engineers in Rust who need a complete, idiomatic toolkit.
• Python teams that want high-performance inference/training through PyO3 bindings with maximum zero-copy.
• Production systems where scikit-learn’s pickle is a security risk and every millisecond counts.

Project Vision

Goal

sciencekit is a Machine Learning library written natively in Rust exposing every algorithm and utility of scikit-learn, rewritten from scratch in an optimized form. It is not a rushed port: each algorithm evolves through mandatory stages — naive implementation, tests, performance tuning (SIMD, parallelism, memory layout) and out-of-core support.

Positioning

Equivalent to scikit-learn in algorithmic coverage, with the performance and safety of Rust. Three target audiences:

Scope

API principles

1. Mandatory builder pattern — every estimator, transformer and pipeline exposes a builder; direct constructors stay private. Every builder accepts execution_mode(SKExecutionMode::...), defaulting to Automatic.
2. Zero-copy inputs — public APIs receive ArrayView, CowArray or sparse views (sprs), never arrays by value.
3. Type-safe pipelines — step compatibility validated at compile time through associated types.
4. Full descriptive names — no abbreviations; the single exception is the project prefix sk/SK. Examples: maximum_number_of_iterations (not max_iter), nearest_neighbors_count (not k).

// Core trait vocabulary — implemented by the sciencekit crates.
pub trait SKEstimator { /* hyperparameters and fit */ }
pub trait SKRegressorPredictor: SKEstimator { /* predict -> Array1<F> */ }
pub trait SKClassifierPredictor: SKEstimator { /* labels -> Array1<i64> */ }
pub trait SKTransformer: SKEstimator { /* transform / fit_transform */ }

pub trait SKDataSource { /* eager: full in-memory access */ }
pub trait SKLazySource { /* streaming: Iterator<Item = Batch> */ }
pub trait SKMappableSource { /* memmap: O(1) random access */ }

Algorithms compose exactly the traits they support: an SKSGDClassifier implements SKClassifierPredictor + SKLazySource; an SKStandardScaler implements only SKTransformer.

Interoperability

• Single internal format: extended Safetensors — JSON header with hyperparameters, training state and recoverability metadata (checkpointing), with partial writes via sharding or header padding.
• ONNX: every estimator implements SKToOnnx; external models load as a generic predictor.
• Data: pluggable sources behind arbitrary I/O traits, with conversions for Polars and DataFusion.
• Observability: tracing oriented to OpenTelemetry, disableable at near-zero cost.

Algorithms

The catalog mirrors the complete scikit-learn taxonomy. Every item enters the roadmap with a mandatory four-stage evolution — naive → tests → performance → out-of-core — and receives a Python binding as soon as it is validated on CPU.

Preprocessing · sciencekit_preprocessing

Imputation · sciencekit_impute

Linear models · sciencekit_linear_model

Nearest neighbors · sciencekit_neighbors

Support vector machines · sciencekit_svm

Trees and ensembles · sciencekit_tree + sciencekit_ensemble

Unsupervised · sciencekit_cluster, sciencekit_decomposition, sciencekit_outlier

Model selection, pipelines and metrics

Architecture

The project is a Cargo workspace split into cohesive sub-crates, optimizing build times and isolating dependencies. Every structural rule below comes from the PRD (§3–§9).

Workspace layout

flowchart TD
    accTitle: The sciencekit Cargo workspace layout
    accDescr: A workspace root holding Cargo.toml and a crates folder. The crates folder contains an umbrella crate, two foundation crates, and one crate per algorithm area, plus a docs folder.
    ROOT[sciencekit/]
    ROOT --> CARGO[Cargo.toml<br/>workspace root]
    ROOT --> CRATES[crates/]
    ROOT --> DOCS[docs/]

    subgraph FOUNDATION[Foundation]
        UMBRELLA[sciencekit/<br/>umbrella crate re-exports]
        COMMON[sciencekit_common/<br/>base traits, types, central errors]
        MATH[sciencekit_math/<br/>linear algebra, BLAS, SIMD, sparse, pairwise]
    end

    subgraph ALGORITHMS[Algorithm crates]
        PREP[sciencekit_preprocessing/<br/>scalers, encoders, polynomial features]
        IMPUTE[sciencekit_impute/<br/>simple and multivariate imputers]
        LINEAR[sciencekit_linear_model/<br/>linear regressions, logistic, SGD]
        NEIGHBORS[sciencekit_neighbors/<br/>KNeighbors, KDTree, BallTree]
        SVM[sciencekit_svm/<br/>SVC, SVR, LinearSVC]
        TREE[sciencekit_tree/<br/>decision trees]
        ENSEMBLE[sciencekit_ensemble/<br/>random forest, bagging, gradient boosting]
        CLUSTER[sciencekit_cluster/<br/>KMeans, DBSCAN, MiniBatchKMeans]
        DECOMP[sciencekit_decomposition/<br/>PCA, TruncatedSVD]
        OUTLIER[sciencekit_outlier/<br/>anomaly detection]
    end

    subgraph SUPPORT[Support crates]
        SELECTION[sciencekit_model_selection/<br/>KFold, train_test_split, GridSearchCV]
        PIPELINE[sciencekit_pipeline/<br/>type-safe SKPipeline, DAGs]
        METRICS[sciencekit_metrics/<br/>evaluation metrics]
        INTEROP[sciencekit_interop/<br/>ONNX, Safetensors, I/O traits]
        GPU[sciencekit_gpu/<br/>SKComputeBackend + OpenCL/CUDA/ROCm]
        PYTHON[sciencekit_python/<br/>PyO3 bindings]
    end

    CRATES --> FOUNDATION
    CRATES --> ALGORITHMS
    CRATES --> SUPPORT

    style ROOT fill:#e0f2fe,stroke:#0284c7,color:#0c4a6e
    style FOUNDATION fill:#f0fdf4,stroke:#16a34a,color:#14532d
    style ALGORITHMS fill:#fef3c7,stroke:#d97706,color:#78350f
    style SUPPORT fill:#ffe4e6,stroke:#e11d48,color:#7f1d1d

Structural rules

Naming convention

No abbreviations anywhere — single exception: the project prefix sk/SK.

// Illustrative — naming convention in action.
let model = SKKMeansClassifierBuilder::new()
    .maximum_iterations(300)      // not max_iter
    .nearest_neighbors_count(5)   // not k
    .build()?;

Core traits (sciencekit_common)

Trait composition is the foundation of everything:

// Core trait vocabulary — implemented by the sciencekit crates.
pub trait SKEstimator { /* hyperparameters and fit */ }
pub trait SKRegressorPredictor: SKEstimator { /* predict -> Array1<F> */ }
pub trait SKClassifierPredictor: SKEstimator { /* labels -> Array1<i64> */ }
pub trait SKTransformer: SKEstimator { /* transform / fit_transform */ }

pub trait SKDataSource { /* eager: full in-memory access */ }
pub trait SKLazySource { /* streaming: Iterator<Item = Batch> */ }
pub trait SKMappableSource { /* memmap: O(1) random access */ }
pub trait SKToOnnx { /* ONNX export */ }
pub trait SKComputeBackend { /* compute device abstraction */ }

Algorithms compose exactly the traits they support: an SKSGDClassifier implements SKClassifierPredictor + SKLazySource; an SKStandardScaler implements only SKTransformer.

Memory management

Higher-order iteration is mandatory: .map(), .map_inplace(), .zip_mut_with(), plus azip!() / par_azip!() and rayon for data parallelism — never manual index loops.

Custom allocators

Expected gains vs glibc on ML workloads: 15–30% allocation throughput. The automatic decision mechanism picks the default per platform and workload.

Out-of-core strategies

Every algorithm declares its access pattern and implements what works:

Sparse arrays are first-class from day one via sprs + ndarray (CSR/CSC/COO) — essential for Lasso, linear SVMs and text classification.

Threads and concurrency

Automatic execution decision

Every builder exposes execution_mode(SKExecutionMode::...). Default: Automatic.

The decision weighs: available memory (sysinfo), dataset size, declared access pattern, CPU cores (available_parallelism) and an optional batch hint. Users can override any parameter in the builder.

Size-aware parallelism

Resolution derives parallelism = min(cores, ceil(unit / grain)), where the work unit is the whole dataset (in-memory), one batch (streaming), or an arbitrary shard (memory-mapped → always the core count). A single-core machine or a sub-grain unit resolves to exactly 1 (sequential, no dispatch overhead). The per-kernel grain (minimum elements per thread to amortize dispatch) is measured by the kernel_calibration criterion bench and recorded beside each kernel with machine + date provenance (see math-kernel).

flowchart LR
    accTitle: How an execution intent becomes a plan and how it is executed
    accDescr: An execution intent and a context are resolved into an execution plan that carries a mode and parallelism. In-memory plans run kernels that honor parallelism; streaming plans use a driver with an I/O thread prefetching into a double buffer and a compute pool updating state; memory-mapped plans run parallel shards.
    Intent[Execution intent] --> Resolve[sk_resolve_execution_plan]
    Ctx[Context: memory / cores / size / grain] --> Resolve
    Resolve --> Plan[SKExecutionPlan: mode + parallelism + batch_size + buffer_depth]
    Plan --> Mode{Mode}
    Mode -- In-memory --> Kernels[Kernels honor parallelism]
    Mode -- Streaming --> Driver[Streaming driver]
    Driver --> IoThread[I/O thread]
    Driver --> Compute[Compute pool]
    IoThread -->|prefetch k+1| Buffer[(Double buffer)]
    Buffer --> Compute
    Compute -->|update batch k| Update[Algorithm update step]
    Update -->|Continue or Stop| Driver
    Mode -- Memory-mapped --> Shards[Parallel shards, parallelism = cores]

For streaming, the driver owns the I/O ∥ CPU overlap: a dedicated I/O thread prefetches the next batch into a bounded double buffer while the compute pool processes the current batch via an algorithm-supplied update(batch, parallelism, state) step. Resolution guards the double-buffer footprint — $2 \times \text{batch} \le \text{available memory}$ — and refuses oversized batch hints before any data is processed.

Optimization layers

1. Algorithmic — correct complexity before micro-optimization.
2. Memory layout — contiguity, cache-friendliness, SoA vs AoS where relevant.
3. Higher-order functions — azip!/par_azip!/zip_mut_with (mandatory).
4. SIMD — LLVM autovectorization → explicit wide on profiled hot paths → std::simd when stable.
5. BLAS — opt-in flag (blas-backend: OpenBLAS default, BLIS alternative, MKL never a direct dependency).
6. GPU — additional backend per separate change via SKComputeBackend (OpenCL → CUDA → ROCm).
7. Allocators — automatic choice, user-configurable.

Feature flags have dual granularity: per capability (parallel, allocator-*, blas-*, gpu-*, telemetry-opentelemetry) and per algorithm group (classification, regression, clustering, …).

Model export and import

• Single internal format: extended Safetensors — JSON header with "sciencekit_format_version", hyperparameters, training state and recoverability metadata (checkpointing).
• Partial writes: in-place tensor updates at known offsets; sharding for large models; header padding for small ones.
• Compression: .safetensors.gz, .safetensors.brotli, .snappy.safetensors.
• ONNX: every estimator implements SKToOnnx; external models load as a generic predictor.
• JSON debug: human-readable serialization through the same central export trait.

Errors use thiserror enums per algorithm with uniform propagation of common errors; observability uses tracing oriented to OpenTelemetry, disableable at near-zero cost.

sciencekit_math — computational kernels and BLAS interface

The sciencekit_math crate is the numeric substrate every downstream algorithm crate builds on (Phase 0.3). It provides the reusable operation kernels, memory-layout helpers, pairwise-distance matrices, sparse products and the BLAS/LAPACK backend abstraction that the algorithm crates consume.

Purpose

• Performance-first kernels honouring PRD §4: higher-order iteration with azip!/par_azip! (never manual index loops), SIMD hot paths, and layout awareness.
• Zero-copy public APIs built on ndarray::ArrayView, ndarray::CowArray and sprs views.
• A uniform linear-algebra abstraction ([SKMathBackend]) with a pure-Rust default (faer) and an opt-in system-BLAS path.

Modules

flowchart LR
    accTitle: Modules of sciencekit_math and their key surfaces
    accDescr: sciencekit_math exposes kernels using azip and par_azip, layout helpers, pairwise distance kernels using a SIMD dot, sparse operations over sprs views, and a backend that defaults to faer with an opt-in blas-backend.
    A[sciencekit_math] --> B[kernels]
    A --> C[layout]
    A --> D[pairwise]
    A --> E[sparse_ops]
    A --> F[backend]
    B --> G["azip! / par_azip! / zip_mut_with"]
    C --> H["sk_memory_layout / sk_force_contiguous"]
    D --> I["wide SIMD dot kernel"]
    E --> J["sprs views"]
    F --> K["SKFaerBackend (pure Rust)"]
    F --> L["blas-backend (opt-in)"]

Key public API

• Kernels — sk_elementwise_transform, sk_binary_combine, sk_axis_sum, sk_scale_in_place.
• Layout — SKMemoryLayout, sk_is_c_contiguous, sk_is_f_contiguous, sk_force_contiguous.
• Distances — sk_squared_euclidean_distance_matrix, sk_euclidean_distance_matrix, sk_manhattan_distance_matrix, sk_cosine_distance_matrix.
• Sparse — sk_csr_dense_product, sk_sparse_product.
• Backend — SKMathBackend (trait), SKFaerBackend, SKMatrixMultiplyBackend, SKNdArrayLinalgBackend (blas-backend only), sk_default_math_backend.

All structs/traits use the SK prefix and free functions the sk_ prefix (PRD §3.4).

Resolved parallelism and grain

Every dense kernel accepts the resolved plan’s parallelism and honours it: 1 runs the sequential azip!/zip_mut_with form with zero dispatch overhead; > 1 dispatches par_azip! (or a per-row/per-chunk parallel reduction) across the rayon compute pool. Parallel results match the sequential reference within floating-point tolerance.

Policy lives in sk_resolve_execution_plan, not in the kernels. Resolution derives parallelism = min(cores, ceil(unit / grain)) from the work-unit size (the whole dataset in memory, one batch when streaming) capped by the core count, and returns exactly 1 for single-core machines or sub-grain units. Each kernel documents a grain — the minimum elements per work unit needed to amortize parallel dispatch — measured by the kernel_calibration criterion bench (cheap/medium/expensive closures × sizes 10¹–10⁷) and recorded beside the kernel with machine + date provenance:

The same grain applies to in-memory and streaming regimes. Re-calibrating on another machine changes only these constants and their provenance; the behaviour contracts above are unaffected.

Usage examples

Elementwise transform and an in-place scale (both take the resolved parallelism):

#![allow(unused)]
fn main() {
use ndarray::{array, Array1};
use sciencekit_math::sk_elementwise_transform;

let input = array![1.0_f64, 2.0, 3.0];
let transformed = sk_elementwise_transform(&input.view(), |x| 2.0 * x + 1.0, 1);
// transformed == [3.0, 5.0, 7.0]
}

Squared-Euclidean distance matrix between two feature sets:

#![allow(unused)]
fn main() {
use ndarray::array;
use sciencekit_math::sk_squared_euclidean_distance_matrix;

let points = array![[0.0_f64, 0.0], [3.0, 4.0]];
let distance = sk_squared_euclidean_distance_matrix(&points.view(), &points.view());
// symmetric, zero diagonal, non-negative
}

A dense matrix product through the pure-Rust default backend:

#![allow(unused)]
fn main() {
use faer::{Mat, MatRef};
use sciencekit_math::{SKFaerBackend, SKMathBackend};

let a = Mat::<f64>::from_fn(2, 2, |i, j| (i * 2 + j) as f64);
let b = Mat::<f64>::identity(2, 2);
let backend = SKFaerBackend::new();
let product = backend.gemm(a.as_ref(), b.as_ref(), 1.0);
}

These examples mirror the crate’s companion *_tests.rs modules, which are the executable source of truth for each function’s behaviour.

The blas-backend feature

By default the crate links no C/Fortran BLAS: heavy algebra runs on the pure-Rust faer backend. Enabling the blas-backend feature swaps dense algebra to ndarray-linalg over the configured blas-src backend (OpenBLAS/BLIS/MKL/Accelerate); faer remains available as a fallback. sk_default_math_backend selects the active backend.

Small fixed-size matrices (nalgebra)

For 2×2/3×3 hot paths (e.g. PCA covariance, affine transforms) where GEMM overhead dominates, nalgebra 0.35 was considered as a small-matrix fallback. Per the wave-plan-foundation decision (ADR #37) it is deferred, not adopted: sciencekit_math exposes no small-matrix hot paths yet, and faer’s small-matrix performance is assumed adequate. The decision is retargeted to W5.4 (PCA), where the 2×2/3×3 covariance path actually lands.

Safety

This is the performance crate: unsafe is confined to the SIMD (wide) and matrix multiply kernels and never appears in the higher-order or layout layers.

Plan & Roadmap

The roadmap follows the scikit-learn mapping faithfully, respecting dependencies between crates. Development methodology is as important as the destination. The wave decomposition below comes from the wave-plan-foundation planning change (OpenSpec) — every downstream change references it as the shared reference.

Iterative evolution per algorithm

No stage may be skipped — each algorithm walks the full path:

Acceptance checklist (every implementation)

• Runs correctly with lots of data and with little data
• Correct under concurrency (automatic + applicable explicit modes)
• Exports model (minimum Safetensors) and produces metrics
• Companion _tests.rs modules covering everything above
• Complete names, no abbreviations (single exception: sk/SK)
• Public items with mandatory prefix; methods unprefixed
• No file over 200 lines without folder-module organization
• Python binding in an associated change; GPU backends in separate changes

Phases

Wave plan and foundational decisions

The wave-plan-foundation planning change decomposes PRD phases 0–7 into waves W0–W7 — roughly 60 downstream changes — so dependencies are respected naturally and no throwaway/placeholder implementations are needed. It also locks six foundational technical decisions that every downstream change references:

Deferred tracker (ADR #37): the nalgebra role for small fixed-size matrices was deferred at W0.3 — not adopted, since sciencekit_math has no small-matrix hot paths yet. Revisit at W5.4 (PCA), adopting nalgebra 0.35 only if faer’s small-matrix performance proves inadequate for 2×2/3×3 covariance.

Cross-cutting (continuous)

How to add an Algorithm

This is the boiler plate. It contains no algorithm — only the skeleton every algorithm in sciencekit must follow, and the discipline that keeps hundred algorithms coherent instead of a drawer of clever one-offs. Every case-study chapter in this part assumes this recipe and refers to it by section number.

Contract level. Everything below is written against the merged master contracts of sciencekit_common and sciencekit_math as of the parallel-kernel-execution wave: SKBuilder/SKBuilderState, SKSupervisedFit/SKUnsupervisedFit/SKFeatureTransformer/SKRegressorPredictor/SKClassifierPredictor, SKDataView/SKTargetView, SKExecutionContext/sk_resolve_execution_plan, sk_run_streaming_driver, SKLazySource/SKMappableSource, the measured-grain kernels and the SKMathBackend abstraction. When a real algorithm crate lands and any signature has moved, this chapter is the one that gets the maintenance edit.

0. The journey at a glance

The whole routine, before the details:

flowchart TD
    accTitle: The algorithm routine from failing test to acceptance
    accDescr: Test-first contract, builder plus estimator, execution planning across three regimes, model export, and PRD section 8.7 acceptance gates.
    classDef buildNode fill:#e0f2fe,stroke:#0284c7,color:#0c4a6e,stroke-width:2px
    classDef planNode fill:#f0fdf4,stroke:#16a34a,color:#14532d,stroke-width:2px
    classDef regNode fill:#fff7ed,stroke:#ea580c,color:#7c2d12,stroke-width:2px
    classDef outsNode fill:#ede9fe,stroke:#7c3aed,color:#3b0764,stroke-width:2px

    T["TDD: failing<br>contract test first"]:::buildNode
    T --> B["Builder + estimator<br>SKBuilderState, validate in build()"]:::buildNode
    B --> Fit["estimator · fit(&self) implements<br>SKUnsupervisedFit |<br>SKSupervisedFit<br>+ SKFeatureTransformer when transformer"]:::buildNode
    Fit --> P{"SKExecution<br>Context<br>grain · access_pattern"}:::planNode
    P --> PLN["sk_resolve_execution_plan<br>(Automatic default)"]:::planNode
    PLN -->|"InProcess* | whole array"| K["measured-grain kernels<br>plan.parallelism threads"]:::regNode
    PLN -->|">RAM, sequential"| Dr["sk_run_streaming_driver<br>State + update() per batch"]:::regNode
    PLN -->|">RAM, RandomAccess"| Ma["SKMappableSource<br>memory-mapped rows"]:::regNode
    K -->     Ms["Model (distinct type)<br>predictor /<br>transform verbs"]:::outsNode
    Dr --> Ms
    Ma --> Ms
    Ms --> Chk["§8.7 acceptance + export + scorers<br>fmt · clippy · test · mdbook build"]:::outsNode

1. What “adding an algorithm” means

You are not writing “a file with a math loop”. You are adding a new public estimator to the catalog, which means six obligations, each checked by a gate:

Skipping any of the six produces the estimators this library exists to make unnecessary.

2. File layout and naming

2.1 Folder modules, always private-state, never index loops

Once any file reaches 200 lines, it becomes a folder module with a standardized shape. Start in that shape directly: implement the algorithm expecting to exceed 200 lines.

crates/sciencekit_preprocessing/src/standard_scaler/
├── mod.rs                     # pure dispatcher: mod declarations + re-exports ONLY
├── builder.rs                 # builder + hyperparameter validation + build()
├── core_implementation.rs     # the estimator struct + fit internals
├── model.rs                   # the fitted object + transform
└── standard_scaler_tests.rs   # companion tests, declared at the END of mod.rs

Three hard rules from the conventions chapter:

1. mod.rs carries no logic. It declares submodules and re-exports public items.
2. Tests live beside the implementation in a companion *_tests.rs — never a global tests/ directory. The tests are later quoted, verbatim, in the Development Explained chapter for your algorithm; they are the executable truth.
3. No manual index loops. Iteration is azip!(), map(), par_azip!().

2.2 Naming rules (PRD §3.4)

• No abbreviations, ever, with the single exception of the project prefix. The builder field is maximum_number_of_iterations: usize, not max_iter; the count of nearest neighbors is nearest_neighbors_count, not k.
• Structs and traits: SK + PascalCase → SKStandardScaler, SKRegressorPredictor.
• Free public functions, public variables, public modules: sk_ + snake_case → sk_resolve_execution_plan, sk_axis_sum.
• Methods get no prefix — standing inside an impl is the privilege that earns it: fn fit(...), fn transform(...), fn predict(...).
• Crates keep their full name: sciencekit_preprocessing, never sk_prep.

The discipline is not pedantry: a dataframe.max_iter reads like a contract to every future reader; an abbreviation makes every editor and greppable name a guessing game.

3. The four contracts your algorithm plugs into

sciencekit does not have one Estimator god-trait. It has small verbs, and your algorithm composes exactly the ones it supports (PRD §3.3).

3.1 Fitting traits — one of three

#![allow(unused)]
fn main() {
// Unsupervised (scalers, kmeans, pca): features only.
pub trait SKUnsupervisedFit<F: SKFloat> {
    type Model;                                  // the FITTED object — a distinct type
    type Error: From<SKError>;
    fn fit<'a, X>(&self, features: X) -> Result<Self::Model, Self::Error>
    where
        X: TryInto<SKDataView<'a, F>, Error = SKError>;
}

// Supervised (regressors, classifiers): features + targets.
pub trait SKSupervisedFit<F: SKFloat> {
    type Model;
    type Error: From<SKError>;
    fn fit<'a, X, T>(&self, features: X, targets: T) -> Result<Self::Model, Self::Error>
    where
        X: TryInto<SKDataView<'a, F>, Error = SKError>,
        T: TryInto<SKTargetView<'a>, Error = SKError>;
}

// Transformers additionally transform, with a statically known output type
// (this associated type is what makes pipeline chaining compile-time checked).
pub trait SKFeatureTransformer<F: SKFloat> {
    type Output;
    type Error: From<SKError>;
    fn transform<'a, X>(&self, features: X) -> Result<Self::Output, Self::Error>
    where
        X: TryInto<SKDataView<'a, F>, Error = SKError>;
}
}

Three details carry the design and are easy to miss in a skim:

• fit(&self, ...) — not &mut self. An estimator is immutable configuration; you may hand the same SKStandardScaler to two threads fitting two datasets concurrently. All mutation lifetime lives in the returned Model.
• Self::Model is a distinct type. Do not reuse the estimator struct as its own model. The model holds fitted values (mean, variance, coefficients…), implements the operation traits, and is what gets serialized for export.
• Everything enters through TryInto<SKDataView> / TryInto<SKTargetView>. The caller passes zero-copy inputs; the fallible conversion seam (Error = SKError, the central taxonomy: ShapeMismatch, UnsupportedRepresentation, …) is where validation lives. Never accept raw Array2 by value.

3.2 The model’s predictive side

#![allow(unused)]
fn main() {
// Implemented on the MODEL, never on the estimator.
pub trait SKRegressorPredictor<F: SKFloat> {
    type Error: From<SKError>;
    fn predict<'a, X>(&self, features: X) -> Result<ndarray::Array1<F>, Self::Error>
    where
        X: TryInto<SKDataView<'a, F>, Error = SKError>;
}

pub trait SKClassifierPredictor<F: SKFloat> {
    type Error: From<SKError>;
    fn predict_labels<'a, X>(&self, features: X) -> Result<ndarray::Array1<i64>, Self::Error>
    where
        X: TryInto<SKDataView<'a, F>, Error = SKError>;
    fn predict_probabilities<'a, X>(&self, features: X) -> Result<ndarray::Array2<F>, Self::Error>
    where
        X: TryInto<SKDataView<'a, F>, Error = SKError>;
}
}

Note predict returns Array1<F> in the model scalar — an f32 model predicts f32 with no f64 detour — while classifier labels are scalar-independent i64 indices (SKTargetView::Continuous is likewise generic over F).

4. The builder

The builder is the only public constructor path. Direct constructors are private; a new estimator validates every hyperparameter exactly once, in build().

#![allow(unused)]
fn main() {
use sciencekit_common::builders::{SKBuilder, SKBuilderState};
use sciencekit_common::builders::sk_validate_hyperparameter;
use sciencekit_common::SKExecutionMode;

pub struct SKStandardScaler {
    plan: SKBuilderState,
    with_mean: bool,
    with_std: bool,
}

impl SKStandardScaler {
    pub fn new() -> SKStandardScalerBuilder {
        SKStandardScalerBuilder::new()
    }
    // ... the estimator also exposes:
    pub fn execution_intent(&self) -> SKExecutionMode { self.plan.execution_intent() }
    pub fn resolve_plan(
        &self, context: &SKExecutionContext,
    ) -> Result<SKExecutionPlan, SKError> { self.plan.resolve_plan(context) }
}

pub struct SKStandardScalerBuilder {
    state: SKBuilderState,     // holds execution_intent, defaults to Automatic
    with_mean: bool,
    with_std: bool,
}

impl SKStandardScalerBuilder {
    pub(crate) fn new() -> Self {
        Self { state: SKBuilderState::new(), with_mean: true, with_std: true }
    }

    pub fn with_mean(mut self, value: bool) -> Self { self.with_mean = value; self }
    pub fn with_std(mut self, value: bool) -> Self { self.with_std = value; self }

    pub fn execution_mode<E>(mut self, mode: E) -> Self
    where E: TryInto<SKExecutionMode> {
        self.state.execution_mode(mode);
        self
    }
}

impl SKBuilder<SKStandardScaler> for SKStandardScalerBuilder {
    fn execution_mode<E>(&mut self, mode: E) -> &mut Self
    where E: TryInto<SKExecutionMode> {
        self.state.execution_mode(mode);
        self
    }
    fn build(self) -> Result<SKStandardScaler, SKError> {
        // validation here — e.g. with_mean=false && with_std=false leaves identity.
        Ok(SKStandardScaler { plan: self.state, with_mean: self.with_mean, with_std: self.with_std })
    }
}
}

The obligations encoded:

• execution_mode(SKExecutionMode::Automatic) is the default resolved intent; the user may override, and the trait’s &mut self form coexists with the by-value chaining form — both funnel into the same SKBuilderState.
• build() never panics; invalid states return SKError::InvalidHyperparameter via sk_validate_hyperparameter(name, valid, reason).

4.1 Why the intent is resolved at operation time, not at build time

A constant that tempts every newcomer: resolve the plan eagerly in build() and store it. Do not. The estimator is reusable: the same fitted model may be asked to transform a 100-row slice in the morning and drive the usual 50 GB job in the afternoon — and available_memory_bytes, the dominant variable of the resolution rule, is free memory right now. Resolution happens per-operation, inside fit/transform/predict, from a fresh SKExecutionContext.

5. Planning: from intent to executable plan

Inside every operation you do the same four lines, and then parametrize everything by plan.parallelism and guard everything by plan.mode:

#![allow(unused)]
fn main() {
use sciencekit_common::execution::{SKExecutionContext, sk_resolve_execution_plan};
use sciencekit_common::observability::{SKOperationAttributes, SKOperationKind};
use sciencekit_math::kernels::SK_AXIS_SUM_GRAIN; // your dominant kernel's grain

fn standardize<F: SKFloat>(
    estimator: &SKStandardScaler,
    features: &ndarray::ArrayView2<F>,
) -> Result<ndarray::Array2<F>, SKError> {
    sk_run_operation(
        SKOperationAttributes {
            operation: SKOperationKind::Fit,
            rows: features.nrows(),
            columns: features.ncols(),
            execution_mode: estimator.execution_intent(),
            backend: SKBackendKind::Faer, // or the actual backend's kind()
        },
        || {
            let mut context = SKExecutionContext::real();
            context.dataset_size_bytes = features.len() as u64 * size_of::<F>() as u64;
            context.scalar_size_bytes = size_of::<F>() as u64;   // defaults are f64-sized; set it for f32 batches
            context.dataset_elements = Some(features.len() as u64);
            context.access_pattern = SKAccessPattern::Sequential;
            context.grain = SK_AXIS_SUM_GRAIN;          // the algorithm declares its grain
            let plan = sk_resolve_execution_plan(estimator.execution_intent(), &context)?;

            // dispatch on plan.mode / use plan.parallelism (see §6, §7)
            // ...
            unimplemented!()
        },
    )
}
}

What sk_resolve_execution_plan does with that (pure, deterministic — so your tests can assert it exactly):

intent = Automatic:
    dataset_size_bytes ≤ available_memory_bytes  →  InProcessSynchronous
    else access_pattern RandomAccess             →  OutOfCoreMemoryMapped
    else Sequential | Iterative                  →  OutOfCoreStreaming

parallelism = min(cores, ceil(unit_elements / grain))
              unit = whole dataset (in-memory) | one batch (streaming, batch_size_hint) | ∞ (mapped)

Two compile-time facts your tests should lock in:

• Automatic is a policy, not a prophecy: it never returns ExecutionModeIncompatible, and small data always lands in-memory (parallelism collapses to 1 when the unit is below one grain — the resolver will not hand out threads for a 5-element array). One caveat: when Automatic resolves to a streaming plan, the buffer guard still applies — an oversized batch_size_hint raises BatchBufferOversize, which is a data-shape error, not an intent error.
• Declaring OutOfCoreStreaming with access_pattern = RandomAccess is a programming error caught before any data moves: SKError::ExecutionModeIncompatible.

5.1 Which access pattern should your algorithm declare?

The pattern is a property of the algorithm math, not a preference. SKStandardScaler is Sequential. SKSGDClassifier is Sequential per pass (it never needs a random row). SKKNeighborsClassifier is the textbook RandomAccess: its “training” is storing and its prediction needs any stored row at any time.

6. Regime 1 — in-memory, parallel kernels

InProcessSynchronous is the regime every algorithm supports: the whole array is resident, and the job is to throw plan.parallelism threads at it. You do not create threads and you do not write rayon code in the algorithm — the kernels do, taking a single parallelism: usize and branching internally:

#![allow(unused)]
fn main() {
use sciencekit_math::kernels::sk_axis_sum;

// axis 0 => per-column sums; parallel path uses per-chunk PARTIAL accumulators
// combined at the end (no thread ever writes to another thread's accumulator).
let per_column_sums = sk_axis_sum(&features.view(), 0, plan.parallelism);
}

The rules to internalize before writing your first kernel call:

• Parallelism is paid for by the grain. Each kernel family carries a measured grain (SK_AXIS_SUM_GRAIN = 100_000 elements for reductions, 1M for elementwise transforms…) — the crossover point where splitting work pays for its own dispatch overhead. You inject your dominant kernel’s grain into the context; the resolver computes parallelism = min(cores, ceil(unit / grain)). The same number governs in-memory and streaming: the grain is a property of the kernel, not of the regime.
• Partials, never shared accumulators. When you compute anything in parallel, each worker accumulates into its own partial vector; partials are combined once at the end. This is why the parallel column sum is deterministic-order-tolerant and race-free — and in the transformer chapter it is why Σ and Σx² partials compose into Welford cleanly.
• Layout is a contract, not an accident. Pass &features.view() into kernels; when a copy may be needed for a blocked/backend path, sk_force_contiguous gives back a CowArray — borrow if already contiguous, own only if strided.

7. Regime 2 — out-of-core streaming, the generic driver

The mental model first. Streaming is not “a smaller fit called repeatedly by the user”. It is the same fit, running over an SKLazySource that yields owned SKDataBatches, with your algorithm’s state updated per batch inside the driver’s prefetch pipeline. You write one closure; the driver owns buffering, an I/O producer thread, ordering, errors and the final-batch flag.

#![allow(unused)]
fn main() {

struct StreamingMomentsState {
    count: u64,
    count_f: f64,          // Welford accumulators — see the transformer chapter
    mean: ndarray::Array1<f64>,
    m2: ndarray::Array1<f64>,
}

fn fit_streaming<F: SKFloat>(
    source: &mut impl SKLazySource<F>,
    plan: &SKExecutionPlan,
) -> Result<StreamingMomentsState, SKError> {
    let mut state = StreamingMomentsState { count: 0, count_f: 0.0,
        mean: ndarray::Array1::zeros(0), m2: ndarray::Array1::zeros(0) };

    sk_run_streaming_driver(source, &mut state, |batch, parallelism, state| {
        // ← THE "partial fit": one batch, one update of the combinable state.
        update_moments(batch.data(), parallelism, state);
        SKStreamDecision::Continue
    }, plan)?;

    Ok(state)
}
}

Mechanics you inherit for free by using the driver (and must never re-implement by hand — that is a whole second book of subtle concurrency):

• a dedicated I/O producer thread reads ahead through a bounded buffer sized by plan.buffer_depth (double buffering = 1), so disk and CPU overlap;
• batches arrive in source order, owned (SKDataBatch owns its Array2<F> precisely so it can be handed across the producer/consumer boundary);
• SKDataBatch::is_final() marks the last batch; position() gives zero-based order;
• a read error mid-stream surfaces as SKError with all earlier batches’ effects preserved in your state — error handling is monotonic: never a partial model without an error, and never an error that discards the caller’s ability to see it;
• SKStreamDecision::Stop halts the pipeline early while keeping the current batch’s effect.

The same state produced by the stream must be the same kind of state the in-memory path produces. For a scaler that is trivial; for LinearRegression it is exactly why Chapter 4 accumulates a Gram matrix per batch; for SGD Chapter 3 of this part never had to re-formulate anything — its update is already the only formulation it has.

Batch sizing. The plan carries plan.batch_size (streaming only) and checks a buffer guard: 2 × batch_elements × scalar_size ≤ available_memory_bytes, else SKError::BatchBufferOversize. Pass batch_size_hint in the context if your format has a natural batch (e.g. rows of a parquet row-group, cache-line multiples of a .npy memmap).

8. Regime 3 — memory-mapped random access (SKMappableSource)

The third regime exists because not all algorithms read data in order. k-NN needs “row arbitrary”; so do k-means over the same dataset used twice, and decision-tree splits over a base too big for RAM.

#![allow(unused)]
fn main() {
pub trait SKMappableSource<F> {
    type Error: From<SKError>;
    fn number_of_rows(&self) -> usize;
    fn row(&self, index: usize) -> Result<ndarray::ArrayView1<'_, F>, Self::Error>;
}
}

Design rule: implement SKMappableSource when the algorithm declares SKAccessPattern::RandomAccess and the data backing it lives on disk in a layout that gives O(1) rows (contiguous row-major, no header tricks). Then sk_resolve_execution_plan(Automatic)’s resolution sends huge datasets to OutOfCoreMemoryMapped where the union of memory and page cache serves the algorithm — no SKDataBatch copying, no prefetch pipeline, just rows handed out from where they already are.

The transform chapter will not need it; the neighbors chapter is built on it.

9. Observability: every public operation under one span

#![allow(unused)]
fn main() {
use sciencekit_common::observability::{sk_run_operation, SKOperationAttributes, SKOperationKind};

pub fn fit<F>(...) -> Result<Model, MyError> {
    sk_run_operation(
        SKOperationAttributes {
            operation: SKOperationKind::Fit,          // Fit · FitTransform · FitPredict · PartialFit · Transform · Predict · Score
            rows, columns,
            execution_mode, backend,
        },
        || { /* the real work */ },
    )
}
}

sk_run_operation opens a span carrying operation/rows/columns/mode/backend, times the closure and, if it errors, records the error — bounded enums only, never free-form strings. The streaming tip: each per-batch update is instrumented with SKOperationKind::PartialFit (k-means over 10² batches, an SGD classifier with 100 → a streaming job becomes a legible sequence of partial steps rather than one opaque fit).

10. The routine, end to end

Assemble everything above into the routine you will now repeat. It is the same routine no matter the algorithm — only the math body changes:

TDD guarantee (mandatory, never skipped):
  1. Failing contract test FIRST         (cargo test → compile-error/fail; commit on the red)
  2. Smallest naive implementation        (passes; commit)
  3. Performance pass                     (grain, layout, partials; same tests still green;
                                           commit)
  4. Streaming/out-of-core pass           (same state, driver wiring; SAME tests, plus batch
                                           tests; commit)

Gates, before any push:
  cargo fmt --all --check && cargo clippy --workspace --all-targets && cargo test --workspace
  mdbook build docs   (documentation change accompanies every code change)

Acceptance, per PRD §8.7 — the four scenarios each estimator must survive:
  ✔ large and small datasets        (the same tests ask both shapes of data from the contract)
  ✔ under concurrency               (Send/Sync by design; parallel kernels; &self fit)
  ✔ model export & reload           (fit→serialize→deserialize→predict must equal)
  ✔ metrics computed                (scorer traits give you score() for free)

Numbers to keep in mind when you pick the second pass (the ones teams constantly trip on):

• Do not skip the performance stage. The naive version exists to make the next test meaningful: the parallel/Welford/laid-out version must produce the same answers to 1e-15 (see the transformer chapter’s cancellation walk-through), and the grain constants in sciencekit_math::kernels::grain are measured, not folklore — reuse them.
• Automatic is the golden default because resolution is pure and total — the plan is a function of (intent, context) with no world effects, which means every algorithm’s dispatch itself is unit-testable without disks or clusters. Keep it that way: never read the clock, the file system, or the thread pool inside a decision path.
• The model type is the export boundary. PRD §8.7 acceptance includes export; design the model as a plain, default-repr struct of fitted values + the metadata needed to rebuild, and serialization is nearly a recording note rather than a new design.

11. Checklist — before you open the PR

• Folder module; mod.rs is semantics-free; tests are companions, not in tests/
• Naming: no abbreviations; SK/sk_ prefixes applied per §3.4; crate name intact
• Builder-only construction; execution_mode default Automatic; validation in build()
• Zero-copy inputs via TryInto<SKDataView>; errors from the central taxonomy
• Estimator is immutable (&self methods), model is distinct type + implements its verbs
• Execution context declares grain and access_pattern; plan resolved per operation
• In-memory path goes through measured-grain kernels with plan.parallelism
• Streaming path (when the algorithm afford sequential batches): same state + driver; memory-mapped path (if RandomAccess): SKMappableSource
• Every entry point wrapped in sk_run_operation with bounded attributes
• Red→green→refactor committed per task; OpenSpec tasks.md - [x] kept in lockstep
• fmt/clippy/test/mdbook build all green locally
• Documentation branch: Development Explained chapter (new) + tutorials walk-through
• Independent review: spec met, then the PR references the ADR issue with closing keywords

You now know everything that is the same across algorithms. The next four chapters are what is different — the math, the numerics, and the traps of each family.

Tutorials

This part is the hands-on textbook for the library. Where Development Explained walks through what was already built and why, this part teaches you how to build the next thing — by hand, with no agent, using only the book as your companion.

Who this is for

You already know Rust (ownership, generics, traits, iterators — where a Rust language concept appears, we reference the Development Explained book instead of re-teaching it). You do not need to know:

• how sciencekit wires its internals (we show every wire),
• the numerical-computing theory behind the algorithms (we teach each one from zero),
• how parallel machines think about work (we teach grain, partials and prefetch pipelines in context, the first time an algorithm needs them).

The one idea to hold

Every estimator in sciencekit follows the same skeleton:

flowchart TD
    classDef buildNode fill:#e0f2fe,stroke:#0284c7,color:#0c4a6e,stroke-width:2px
    classDef fitNode fill:#f0fdf4,stroke:#16a34a,color:#14532d,stroke-width:2px
    classDef regNode fill:#fff7ed,stroke:#ea580c,color:#7c2d12,stroke-width:2px
    classDef outNode fill:#ede9fe,stroke:#7c3aed,color:#3b0764,stroke-width:2px

    B["Builder::new()<br>hyperparameters · execution_mode(…)<br>(default: Automatic)"]:::buildNode
    B -->|"build()"| E["Estimator<br>(stateless config, &self everywhere)"]:::buildNode
    E -->|"fit(features[, targets])"| M["Model<br>(fitted, shareable, exportable object)"]:::fitNode
    M --> P{"one resolved plan<br>decides the regime"}:::regNode
    P -->|"InProcessSynchronous"| A["parallel kernels, whole array,<br>plan.parallelism threads"]:::regNode
    P -->|"OutOfCoreStreaming"| S["sk_run_streaming_driver:<br>owned SKDataBatch,<br>update(batch, parallelism, &mut State)"]:::regNode
    P -->|"OutOfCoreMemoryMapped"| Mm["SKMappableSource:<br>O(1) random rows,<br>no materialization"]:::regNode
    A --> O["transform / predict verbs"]:::outNode
    S --> O
    Mm --> O

The key that makes all of this cohere: the algorithm is first designed as a function of a combinable state (Σ sums, Welford moments, a gradient accumulator, a set of stored neighbor vectors). Once the state exists, the in-memory path is “the streaming path with one batch that happens to hold everything” — the two regimes share the same math, and the user never picks a code path, only a SKExecutionMode — usually the default Automatic, which resolves it from dataset size, free memory and access pattern.

Chapter map

Each case-study chapter follows the same internal contract:

1. Story & decisions — what problem, which roads taken and not taken.
2. Theory from zero — numerics with sources; pros/cons; where the obvious solution quietly fails.
3. A state that combines — the mathematical core both regimes share.
4. Full reference code — real, final-form files (folder module, builder, estimator, model, tests).
5. Both regimes side by side — Mermaid data flow and when to prefer each.
6. Acceptance walk-through — PRD §8.7 mapped to the tests you just wrote.
7. Study sources.

Prerequisites check

Before starting, make sure you can run the CI gates locally — they are part of every task:

$ cargo fmt --all --check
$ cargo clippy --workspace --all-targets
$ cargo test --workspace

and build the book you are reading:

$ mdbook build docs

Everything else — grain constants, error taxonomy, planning rules, backend selection — is introduced where the algorithms first need it.

Adding a Transformer — SKStandardScaler and SKRobustScaler

Pre-read: How to add an Algorithm — this chapter quotes its skeleton (folder layout §2, builder §4, planning §5, streaming driver §7) without repeating it.

A transformer is the smallest complete algorithm: it has fit, it has transform, it has no targets, and it has a numerical half-life — one where the naive solution is still the answer printed in most textbooks and, precisely for that reason, is the wrong answer to ship first.

Everything in one place: after this chapter you can write a transformer end-to-end by hand.

flowchart LR
    classDef planNode fill:#e0f2fe,stroke:#0284c7,color:#0c4a6e,stroke-width:2px
    classDef computeNode fill:#f0fdf4,stroke:#16a34a,color:#14532d,stroke-width:2px
    classDef streamNode fill:#fff7ed,stroke:#ea580c,color:#7c2d12,stroke-width:2px
    classDef modelNode fill:#ede9fe,stroke:#7c3aed,color:#3b0764,stroke-width:2px

    B["SKStandardScalerBuilder<br>.execution_mode(...)"]:::planNode --> E["SKStandardScaler<br>(immutable estimator)"]:::planNode
    E -->|fit x| R{"sk_resolve_execution_plan<br>(Automatic)"}:::planNode
    R -->|fits in RAM| M["in-memory regime<br>sk_axis_sum partials<br>parallelism threads"]:::computeNode
    R -->|larger than RAM| S["streaming regime<br>sk_run_streaming_driver<br>Welford update per batch"]:::streamNode
    M --> St["SKStandardScalerModel<br>μ, σ² (fitted)"]:::modelNode
    S --> St
    St -->|transform| T["sk_scaled features"]:::modelNode

1. Story & decisions

The problem. Your features live on different scales — income in the tens of thousands, age in tens. Any distance- or gradient-based model becomes a physics experiment where one feature’s magnitude is the whole of gravity. Standardization removes feature scale as a variable: subtract the mean, divide by the standard deviation. Robust standardization does the same with order statistics (median, inter-quartile range) instead of moments, because a single bad sensor reading can drag a mean half a column away.

The contract shapes the math

SKFeatureTransformer::Output must be statically known — that is the compile-time pipeline hook. For both scalers the output is a materialized owned Array2<F> on transform (zero-copy input goes in, owned result comes out; in-place variants are an explicit future addition, see §5).

2. Theory from zero — moments, and why the obvious solution fails

2.1 The mean, and what “sum” quietly means

The mean of column j is μ_j = (1/N) Σᵢ xᵢⱼ. The variance is σ²_j = (1/N) Σᵢ (xᵢⱼ − μ_j)². Both look like one-liners. They are one-liners in exact arithmetic — and the whole content of this chapter lives in what floating-point arithmetic does to the naive forms.

Parallelism. Sums are associatively commutative in exact math and not quite in floating point — which is why column sums are computed as per-chunk partials and combined once (sk_axis_sum, axis 0). Every parallel worker owns a partial; the combine step happens once. Mutating a shared accumulator from several threads is a data race — and, even if you patched the race with an atomic, you would get a non-deterministic summation order that changes run to run (deviating between small and large data, the exact thing the acceptance §8.7 “little and lots” tests want to distinguish: with little data, the deployment-time parallelism is 1; with lots, the parallel path must agree to machine precision).

2.2 Catastrophic cancellation — the three-line formula that lies

The naive algebraic identity for variance is:

σ² = Σx²/N − μ²

Pro: one accumulation, reads the array once, easy to parallelize (Σ and Σx² partials in the same pass). Con: for data whose mean is huge compared to its spread — a very common shape in financial data, sensor data, any data conditioned by adding a constant — you subtract two numbers that are nearly equal, each with ~16 significant decimal digits of a 64-bit float, and the difference carries only the digits the two big numbers did not have in common.

A classic reminder (need to show scale, not precision):

σ² = Σx²/N − μ²    where each summand carries 53 bits
|Σx²/N − μ²|  ≪  Σx²/N          ⇒ information DESTROYED has exactly the parts you needed

with a proviso: when the spread is ~1e-8 or less of the mean, the naive identity is not “slightly imprecise” — it returns negative numbers that square to imaginary standard deviations (a real bug class reported repeatedly in scikit-learn; see the “welford-negative-variance” folklore and the RunningStats discussions).

2.3 Welford’s online algorithm — the combining trick

The fix (Welford 1962; the multi-participant combiner is Chan, Golub & LeVeque 1983, “Algorithms for computing the sample variance: analysis and recommendations”, The American Statistician) re-represents the state as (count, mean, M2) rather than (n, Σx, Σx²):

W(x; count, mean, M2):
    count     = count + 1
    delta     = x − mean
    mean     = mean + delta / count          // new mean = old mean + (x − old mean)/n
    M2       = M2 + delta·(x − mean)          // M2 = Σᵢ (xᵢ − current_mean)²

output: variance = M2 / count

Every intermediate is a small residual — the possibly-conditioned subtraction is (x − mean), not (big − big); the cancellation is structured, not pathological.

And when the participants are batches (streaming), the form composes: two partial states (count_a, mean_a, M2_a), (count_b, mean_b, M2_b) combine with

delta      = mean_b − mean_a
count_c    = count_a + count_b
mean_c    = mean_a + delta · count_b / count_c
M2_c      = M2_a + M2_b + delta² · (count_a · count_b / count_c)

That is the whole algorithm: for the streaming path the “update per batch” is exactly the above with batch=(count_b, mean_b, M2_b) — each batch collapsed to its OWN Welford state (using the parallel-partials trick within the batch), then combined with the running state. One formulation, two regimes; that is the design point the anchor chapter promised.

Sources to learn this properly:

• Donald Knuth, TAOCP vol. 2, §4.2.2 — where most people first meet Welford (as exercise: hardly the bibliographic patient)
• Chan, Golub, LeVeque (1983), Am. Stat. — the analysis that says why it is stable
• Welford’s original note: Technometrics 4:3 (1962), “A note on a method for calculating corrected sums of squares and products”

3. The combinable state

#![allow(unused)]
fn main() {
// crates/sciencekit_preprocessing/src/standard_scaler/core_implementation.rs
use sciencekit_common::SKFloat;
use sciencekit_math::kernels::sk_axis_sum;

/// A combinatorial object: in-memory partials and streaming batches
/// both arrive as these, and both reduce by the same merge.
#[derive(Debug, Clone, Default)]
pub struct SKMoments {
    pub count: u64,
    pub mean: ndarray::Array1<f64>,
    pub m2: ndarray::Array1<f64>,   // Σ (x − μ)² accumulated around the current mean
}

impl SKMoments {
    /// Pull one *chunk of rows* into the state. `parallelism` comes from the plan:
    /// the chunk's column sums are the parallel partials of `sk_axis_sum`.
    pub fn absorb_chunks<F: SKFloat>(
        &mut self,
        chunks: &ndarray::Array2<F>,
        parallelism: usize,
    ) {
        // ONE pass, no Σx² anywhere
        let chunk_rows = chunks.nrows();
        for row in chunks.rows() {
            let row_f64: ndarray::Array1<f64> = row
                .iter()
                .map(|v| v.to_f64().unwrap_or(0.0))
                .collect();
            // Welford single-element update against the CURRENT mean
            self.count += 1;
            // (Chosen for the parallel column-mean first — the streaming variant below
            //  does the whole-batch combine.)
        }
        // In practice, this inner piece merges the batch-level Welford done by
        // `absorb_batch` (below) so both regimes collapse to one combiner.
    }
}
}

The point of showing the will-not-compile version first is pedagogical: the naive habit is to reach for “two parallel sums”. Do not design the first version that way. The one combinable object you write in practice is:

#![allow(unused)]
fn main() {
impl SKMoments {
    /// parallel-friendly: absorb a whole batch by aggregating the batch's own Welford state
    /// first (so parallel partials cannot race), then COMBINING into the running state.
    pub fn absorb_batch<F: SKFloat>(
        &mut self,
        batch: &ndarray::ArrayView2<F>,
        parallelism: usize,
    ) {
        let batch_size = batch.nrows() as u64;
        if batch_size == 0 { return; }
        let per_batch = SKMoments::welford_of_batch(batch, parallelism);
        self.combine(per_batch);
    }

    fn combine(&mut self, other: SKMoments) {
        // The running state starts empty (`SKMoments::default()` carries
        // zero-length vectors): adopting the first partial outright avoids
        // subtracting mismatched shapes on the very first merge.
        if self.count == 0 { *self = other; return; }
        if other.count == 0 { return; }
        let total = self.count + other.count;
        let delta_w: ndarray::Array1<f64> = &other.mean - &self.mean;
        let running_weight = self.count as f64;
        let borrowed_weight = other.count as f64;
        let total_f = total as f64;
        self.mean = &self.mean + &delta_w.mapv(|d| d * borrowed_weight / total_f);
        let other_m2_delta: ndarray::Array1<f64> = delta_w.mapv(|d| d * d);
        self.m2 = &self.m2 + &other.m2
            + other_m2_delta.mapv(|d| d * (running_weight * borrowed_weight / total_f));
        self.count = total;
    }

    fn welford_of_batch<F: SKFloat>(
        batch: &ndarray::ArrayView2<F>,
        parallelism: usize,
    ) -> SKMoments {
        // The batch itself is reduced via per-chunk partial Welford states
        // (sk_axis_sum-shaped: chunk -> own state -> combine), so parallelism
        // never writes the same accumulator twice.
        let (rows, cols) = batch.dim();
        if rows == 0 { return SKMoments::default(); }
        let chunk = rows.div_ceil(parallelism.max(1));
        let partials: Vec<SKMoments> = batch
            .axis_chunks_iter(ndarray::Axis(0), chunk)
            .into_par_iter()
            .map(|piece| {
                let mut state = SKMoments {
                    count: 0,
                    mean: ndarray::Array1::zeros(cols),
                    m2: ndarray::Array1::zeros(cols),
                };
                for row in piece.rows() {
                    // Single-element Welford update, one row at a time: each
                    // row is ONE observation, so the shared observation counter
                    // advances once per row while every column updates its own
                    // (mean, m2) pair with that same counter value.
                    state.count += 1;
                    let count_as_float = state.count as f64;
                    for (mean, m2, value) in itertools::izip!(&mut state.mean, &mut state.m2, row.iter()) {
                        let next = value.to_f64().unwrap_or(0.0);
                        let delta = next - *mean;
                        *mean += delta / count_as_float;
                        *m2 += delta * (next - *mean);
                    }
                    // Parallel safety comes from the outer split: each rayon
                    // worker owns its private `state`; partials merge once via
                    // `combine` below, so no two threads share an accumulator.
                }
                state
            })
            .collect();
        let mut whole = SKMoments::default();
        for partial in partials { whole.combine(partial); }
        whole
    }
}
}

One caveat that keeps this honest: the exact Welford row-deceleration must be generic over both fact layouts, f32 and f64; state carries f64 internally always — variance accumulated in f32 Σx²-adjacent forms collapses exactly in the case Welford exists to prevent. This is a deliberate micro-decision: mean/variance at f32 costs you the 1e-7 precision the rest of the pipeline relies on. The features stay F; the statistics are f64 unconditionally.

4. Both regimes — the rest of the crown code

#![allow(unused)]
fn main() {
// in-memory path: one batch is the whole array; the driver is unnecessary.
impl SKUnsupervisedFit<F> for SKStandardScaler {
    type Model = SKStandardScalerModel;
    type Error = SKError;

    fn fit<'a, X>(&self, features: X) -> Result<Self::Model, Self::Error>
    where X: TryInto<SKDataView<'a, F>, Error = SKError> {
        let view = features.try_into()?.as_dense()?; // reject sparse up front
        sk_run_operation(
            SKOperationAttributes {
                operation: SKOperationKind::Fit,
                rows: view.nrows(),
                columns: view.ncols(),
                execution_mode: self.execution_intent(),
                backend: SKBackendKind::Faer,
            },
            || {
                let mut context = SKExecutionContext::real();
                context.dataset_size_bytes = view.len() as u64 * size_of::<F>() as u64;
                context.dataset_elements = Some(view.len() as u64);
                context.access_pattern = SKAccessPattern::Sequential;
                context.grain = SK_AXIS_SUM_GRAIN;
                let plan = sk_resolve_execution_plan(self.execution_intent(), &context)?;

                let mut state = SKMoments::default();
                state.absorb_batch(&view, plan.parallelism);
                Ok(SKStandardScalerModel::from_moments(state, self))
            },
        )
    }
}
}

and the streaming twin — note that only the plumbing changes:

#![allow(unused)]
fn main() {
impl SKStandardScaler {
    pub fn fit_streaming_source<S: SKLazySource<F> + Send>(
        &self,
        source: &mut S,
        plan: &SKExecutionPlan,
    ) -> Result<SKStandardScalerModel, SKError> {
        let mut state = SKMoments::default();
        sk_run_streaming_driver(source, &mut state, |batch, parallelism, state| {
            sk_run_operation(
                SKOperationAttributes {
                    operation: SKOperationKind::PartialFit,
                    rows: batch.data().nrows(), columns: batch.data().ncols(),
                    execution_mode: self.execution_intent(),
                    backend: SKBackendKind::Faer,
                },
                || { state.absorb_batch(&batch.data(), parallelism); }
            ).map_err(|_| SKError::Io(std::io::Error::new(std::io::ErrorKind::Other, "batch update failed")))?;
            SKStreamDecision::Continue
        }, plan)?;
        Ok(SKStandardScalerModel::from_moments(state, self))
    }
}
}

The model — spelled out once, because it is the object the “smallest complete algorithm” must finish with (per-column fitted statistics as f64, plus the construction flags, which transform needs to know what the caller asked for):

#![allow(unused)]
fn main() {
pub struct SKStandardScalerModel {
    means: ndarray::Array1<f64>,         // per column; 0.0 where with_mean == false
    scales: ndarray::Array1<f64>,        // 1/σ per column; 1.0 where with_std is false
}

impl SKStandardScalerModel {
    pub(crate) fn from_moments(
        moments: SKMoments,
        options: &SKStandardScaler,      // with_mean / with_std flags
    ) -> Self {
        // Population variance (÷N), matching scikit-learn's default.
        assert!(moments.count > 0, "cannot fit on zero rows — reject empty input in fit()");
        let variance = &moments.m2 / moments.count as f64;
        // `scales` stores the MULTIPLICATIVE factor 1/σ: transform centers
        // then multiplies by it (never divides — division by a tiny σ
        // amplifies rounding; multiplication by a precomputed reciprocal is
        // both faster and the convention `sk_scale_in_place` expects).
        // Columns the caller disabled resolve to identities here, so
        // `columns()` below can hand `transform` ready-to-use triples.
        let scales = variance.mapv(|v| if options.with_std { 1.0 / v.sqrt() } else { 1.0 });
        let means = if options.with_mean { moments.mean } else {
            ndarray::Array1::zeros(moments.mean.len())
        };
        Self { means, scales }
    }
    /// (mean, scale) per column; the flags the caller disabled resolve to identities
    pub fn columns(&self) -> Vec<(usize, f64, f64)> {
        /* zip columns with (0.0, 1.0) wherever the construction flags disabled the step */
    }
}
}

transform (§ below) is a pure flat map — in the model, not the estimator:

#![allow(unused)]
fn main() {
impl<F: SKFloat> SKFeatureTransformer<F> for SKStandardScalerModel {
    type Output = ndarray::Array2<F>;
    type Error = SKError;

    fn transform<'a, X>(&self, features: X) -> Result<Self::Output, Self::Error>
    where X: TryInto<SKDataView<'a, F>, Error = SKError> {
        let view = features.try_into()?.as_dense()?;
        let mut out = view.to_owned();
        for (column_index, mean, scale) in self.columns().into_iter() {
            // `scale` is the precomputed reciprocal 1/σ (or 1.0 where the
            // flags disable scaling): center, then MULTIPLY — never divide.
            let mut column = out.column_mut(column_index);
            for value in column.iter_mut() {
                let centered = value.to_f64().unwrap_or(0.0) - mean;
                *value = F::from_f64(centered * scale).unwrap_or(F::zero());
            }
        }
        Ok(out)
    }
}
}

(prose in real code would use sciencekit_math::kernels::sk_scale_in_place per column with the plan parallelism for the design; here it is condensed to fit prose format — the chapter’s walk-through in the tests applies the actual kernel.)

5. SKRobustScaler — what breaks when “mean” becomes “median”

The same skeleton with one state change and one honest confession.

The same. The estimator/model split, the builder, the zero-copy entry, the plan wiring, SKFeatureTransformer on the model, the streaming driver. Nothing changes structurally.

The change. The state is not moments; it is three quantiles per column — median (50%) and the first/third quartile edges bounding the IQR:

scaled = (x − median) / IQR

The honest conversation. An exact quantile is an order statistic: it needs all the values sorted. Streaming resistant numbers to a floating-point pipeline in one pass is a research-grade set of algorithms:

So the practical rule this chapter hands you: a robust transformer either fits its data or it says so. SKRobustScaler keeps the exact-window implementation and, in streaming mode, requires declaring a maximum window size as a hyperparameter (quantile_window_rows), rejecting construction otherwise — rather than degrading quietly to P². That is the streaming trade-off stated, not hidden: you must either touch the data twice (exact) or pick an estimator and document its error.

Sources for quantiles in a stream:

• Jain & Chlamtac, “The P² algorithm for dynamic calculation of quantiles and histograms without storing observations”, CACM 28(10), 1985 — the famous O(1)-memory one-pass quantile estimator
• Greenwald & Khanna, “Space-efficient online computation of quantile summaries”, SIGMOD 2001 — the ε-accurate predecessor of modern sketches
• scikit-learn, sklearn.preprocessing.RobustScaler documentation — the behavior contract this chapter mirrors (quantile range default, centering/scaling flags)

6. The tests — your executable source of truth

The chapter pattern from the anchor: pick the tests that document the numerics, then the contract ones. Real fixtures as ndarray arrays:

#![allow(unused)]
fn main() {
// standard_scaler_tests.rs (companion; lives beside core_implementation.rs)
#[test]
fn sk_mean_variance_matches_closed_form() {
    // little dataset, f64 — in-memory parallel ps=1 
    let x = array![[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]];
    let scaler = SKStandardScaler::new().build().unwrap();
    let model = scaler.fit(x.view()).unwrap();
    let (mu0, mu1) = (model.mean()[0], model.mean()[1]);
    assert!((mu0 - 3.0).abs() < 1e-12);
    // population variance ε N=3 form, 1 ULP tolerance against the textbook values
    assert!((model.variance()[0] - 8.0 / 3.0).abs() < 1e-14);
    assert!((model.variance()[1] - 8.0 / 3.0).abs() < 1e-14);
}

#[test]
fn welford_matches_two_pass_within_machine_epsilon() {
    // adversarial shape: huge mean, tiny spread — the naive Σx² identity FAILS here
    let base = 1.0e9_f64;
    let x: Vec<f64> = (0..1000).map(|i| base + (i as f64) * 1e-3).collect();
    let mu = x.iter().sum::<f64>() / x.len() as f64;
    let naive = x.iter().map(|v| v * v).sum::<f64>() / x.len() as f64 - mu * mu;
    // naive VIOLATES this assertion (negative variance): document it, never rely on it
    assert!(naive < 0.0);                          // the fixed trap, asserted on purpose 
    let scaler = SKStandardScaler::new().build().unwrap();
    let m = scaler.fit(ndarray::ArrayView2::from_shape((x.len(), 1), &x).unwrap()).unwrap();
    let two_pass_var = x.iter().map(|v| (v - mu) * (v - mu)).sum::<f64>() / 1000.0;
    assert!((m.variance()[0] - two_pass_var).abs() < two_pass_var * 1e-13);
    assert!(m.variance()[0] > 0.0, "welford never returns negative variance");
}

#[test]
fn streaming_matches_in_memory_on_thin_batches() { /* identical fixtures through
    fit_streaming, batched 3-rows-at-a-time; answers agree to 1e-14; the "is_final"
    flag tested by construction (driver guarantees terminal batch delivery) */ }

#[test]
fn concurrency_several_thread_model_transform_output_safe() { /* model is Sync: cloned references
    transformed on rayon, states all agree; represents acceptance "under concurrency" */ }
}

The tests are labeled with their pedagogical role: the first makes sure the math is right once, the second documents why naive variance is not correct here (the genuinely non-obvious transformation decision of this chapter), the third locks the dual-regime promise the PRD asks of every algorithm, and the fourth the concurrency acceptance.

7. Acceptance walk-through & study sources

The PRD §8.7 four-item acceptance, mapped:

Study sources for this chapter:

• Chan, Golub, LeVeque (1983) — the combinating variance analysis (see §2.3’s cites)
• Knuth, The Art of Computer Programming II, §4.2.2 — where Σx² collapse is discussed
• practical style blogs: John D. Cook, “Accurately computing running variance” (the modern favorite walk-through of the algorithm written for engineers)
• scikit-learn: sklearn.preprocessing.StandardScaler and RobustScaler docs and their fit source — the semantics you are reimplementing and the exactness conventions you inherit
• Press et al, Numerical Recipes §14.1 — moments and their pitfalls, the “interview answer” reference

Next chapter turns to the first supervised estimator — one where the whole question changes from “moment of a column” to “which vector of coefficients explains a target”.

Adding a Linear Model — SKLinearRegression

Pre-read: How to add an Algorithm (§5 planning, §7 streaming) and Adding a Transformer (the moments state that combines — keep it in mind; the state here is its first supervised cousin).

Linear regression looks like the easiest algorithm in the catalog and is the best teacher of the two most important numerical-preserving decisions in the whole library: which decomposition carries the solve, and what a streaming formulation of a dense linear system means — including the trade-off that makes it not just a “divided fit”.

1. Story & decisions

Mermaid, both regimes at a glance:

flowchart TB
    classDef planNode fill:#e0f2fe,stroke:#0284c7,color:#0c4a6e,stroke-width:2px
    classDef jacNode fill:#f0fdf4,stroke:#16a34a,color:#14532d,stroke-width:2px
    classDef streamNode fill:#fff7ed,stroke:#ea580c,color:#7c2d12,stroke-width:2px
    classDef modelNode fill:#ede9fe,stroke:#7c3aed,color:#3b0764,stroke-width:2px

    Q{"sk_resolve_execution_plan(Automatic)"}:::planNode
    Q -->|"fits in RAM"| X["A = [1‖X]<br>SVD → pinv(A)·y<br>lstsq via faer backend<br>rank + singular_values surfaced"]:::jacNode
    Q -->|"larger than RAM"| GR["per batch:<br>G := G + BᵀB<br>c := c + Bᵀyᵦ<br>row_count := row_count + batch.nrows()"]:::streamNode
    GR --> END["solve(G, c) at stream end<br>⚠ doubles condition number — stated trade-off"]:::jacNode
    X -->|predict| St["SKLinearRegressionModel<br>coefficients, intercept, rank"]:::modelNode

2. Theory from zero — least squares, and why one line hides a lie

The problem. We have rows xᵢ ∈ ℝᵖ and targets yᵢ ∈ ℝ. We want a vector of coefficients β and a scalar intercept b₀ so that xᵢ·β + b₀ ≈ yᵢ. “Best” is least-squares: minimize Σᵢ (yᵢ − xᵢ·β − b₀)². With the augmented design matrix A = [1 | X] (ones in the first column), that is one matrix equation: minimize ‖Aβ − y‖₂.

Three roads to β, and the difference between them is not style — it is numerics:

Road 1 — Normal equations. Set the gradient of ‖Aβ − y‖₂² to zero: AᵀAβ = Aᵀy. One triangular solve afterwards. The seductive part: AᵀA is smaller than A when rows ≫ columns (a p×p object instead of n×p), which made it the pedagogical favorite for decades. The catch is measurable and unforgiving:

κ(AᵀA) = κ(A)² — the condition number squares.

Condition number is the ratio by which a linear solve can amplify input rounding; counting doubles means doubling the exponent of the loss: if A has condition 10⁷, AᵀA has 10¹⁴, and a 64-bit float has 16 digits — with AᵀA the solve is allowed about two digits of trust at exactly the kind of gritty, correlated-feature datasets scikit-learn gets asked to process daily. And near-multicollinearity is not exotic: it is any pair of features that happen to co-move (age and years-of-experience; two sensors reading the same room).

Road 2 — QR. Orthogonalize A (A = QR), then β = R⁻¹Qᵀy. κ(QR) = κ(A) — the condition number is not squared, because the orthogonal factor is perfectly conditioned by construction. Stable, O(np²), the classic library workhorse.

Road 3 — SVD (the mvp of the “if it is important enough to be in the catalog, it must be stable enough to be in the intro” school). A = UΣVᵀ, so β = Σᵢ σᵢ≠0 (uᵢ·y / σᵢ) vᵢ — the pseudoinverse solution. It costs more than QR and does something QR does not: it shows you the rank. Singular values near zero declare the exact features that make the problem degenerate, and the pseudoinverse path (pinv, rcond = ε·max(m,n) in our backend) is exactly how a least-squares solve degrades gracefully when A is not full-rank, returning (i) the minimum-norm solution and (ii) the information that it did so.

The library’s decision: go through SVD — exposed as the backend’s lstsq, which returns solution, rank, singular_values and residual sums — with QR available for the day a profile demands a cheaper decomposition. One product line (feature blas-backend) selects ndarray-linalg vs faer; the algorithm never sees either (§5 anchor — backends are injectable, dispatch-level truths).

Sources to learn the numerics:

• Golub & Van Loan, Matrix Computations — chapters on least squares, condition numbers, and pseudoinverse: the canonical reference
• Trefethen & Bau, Numerical Linear Algebra — lectures 4–5 and 11–19 build SVD and least squares in a way a working engineer reads in an afternoon
• Higham, Accuracy and Stability of Numerical Algorithms — if you want to prove the claims about normal equations; the chapter’s warnings are his §20

2.1 The streaming road

Streaming sits on the same skeleton — but the naive intuition here doubles down on the same hole. The tempting “solve least squares per batch and average the β” is almost never right: a weighted average hides how the covered rows distributed, and per-batch ranks do not compose.

What composes is the raw material of the normal equations — and that is the honest finding: the streaming regions re-shape the choice to Gram:

G := AᵀA        (p × p, symmetric)
c := Aᵀy        (p, right-hand side)

Each update adds a rank-|B| contribution: for each batch B (rows × p),

G := G + BᵀB
c := c + Bᵀy_batch

At the end, one solve of Gβ = c obtains the coefficients. Exactly one extra fact dumps the price tag on the table: we computed AᵀA, so we are square-conditioning again. Thus the streaming trade-off reads, honestly:

This is also why the streaming chapter of the optimizer family (next chapter) prefers SGD on day one: SGD never forms AᵀA, and its engine turns streaming from a constraint into the algorithm’s natural home.

3. The state — Gram, and what to carry besides

#![allow(unused)]
fn main() {
// crates/sciencekit_linear_model/src/linear_regression/core_implementation.rs
#[derive(Debug, Clone, Default)]
pub struct SKGramState<const WITH_INTERCEPT: bool> {
    row_count: u64,
    gram: ndarray::Array2<f64>,  // AᵀA, (p + intercept) × (p + intercept)
    rhs: ndarray::Array1<f64>,   // Aᵀy,  (p + intercept)
}
}

Notes that carry design weight:

• statistics once again accumulate in f64 unconditionally, like the moments state — when you accumulate BᵀB products, f32×f32→f32 products cross 3 magnitudes of drift per row in adversarial inputs, and the covariates’ first digits are exactly what matters;
• const WITH_INTERCEPT: bool carries through; a 1 turns into a prepended column that the gram absorb adds by chain rule, no separate pass needs reading the data;
• row_count is needed for the residual computation only, not the solve.

4. The crown code, both regimes

#![allow(unused)]
fn main() {
impl<F: SKFloat> SKSupervisedFit<F> for SKLinearRegression {
    type Model = SKLinearRegressionModel;
    type Error = SKError;

    fn fit<'a, X, T>(&self, features: X, targets: T) -> Result<Self::Model, SKError>
    where
        X: TryInto<SKDataView<'a, F>, Error = SKError>,
        T: TryInto<SKTargetView<'a>, Error = SKError>,
    {
        let view = features.try_into()?.as_dense()?;
        let target = targets.try_into()?.as_continuous()?.into_owned();

        sk_run_operation(
            SKOperationAttributes {
                operation: SKOperationKind::Fit,
                rows: view.nrows(),
                columns: view.ncols(),
                execution_mode: self.execution_intent(),
                backend: SKBackendKind::Faer,
            },
            || {
                // 1 · PLAN — same four lines as every algorithm
                let mut context = SKExecutionContext::real();
                context.dataset_size_bytes = view.len() as u64 * size_of::<F>() as u64;
                context.dataset_elements = Some(view.len() as u64);
                context.access_pattern = SKAccessPattern::Sequential;
                context.grain = SK_BINARY_COMBINE_GRAIN; // closeness-of-fit of the batch
                                                          // kernel; the GEMM grain is not yet
                                                          // calibrated — see the note below
                let plan = sk_resolve_execution_plan(self.execution_intent(), &context)?;

                // 2 · SOLVE — through the backend; vanilla-f64 design matrix
                let backend: Box<dyn SKMathBackend<F>> = Box::new(SKFaerBackend::new());
                let (design, target_f) = self.build_design(view, target);
                let least_squares = backend.lstsq(
                    design.view(),              // A, n × p (+1)
target.view().insert_axis(Axis(1)),  // y as n × 1
                )?;

                // 3 · SURFACE — rank and singulars are information, not noise
                SKLinearRegressionModel::from_least_squares(least_squares, self)
            },
        )
    }
}
}

and the streaming twin, which really is “the same data pages, a different sink”:

#![allow(unused)]
fn main() {
impl SKLinearRegression {
    pub fn fit_streaming_source<S, F>(&self, source: &mut S, plan: &SKExecutionPlan)
    -> Result<SKLinearRegressionModel, SKError>
    where
        F: SKFloat + Send,
        S: SKLazySource<F> + Send,
    {
        let mut state = SKGram::zeros(design_width);   // AᵀA & Aᵀy carry the regression
        sk_run_streaming_driver(source, &mut state, |batch, parallelism, state| {
            sk_run_operation(
                SKOperationAttributes {
                    operation: SKOperationKind::PartialFit,
                    rows: batch.data().nrows(), columns: batch.data().ncols(),
                    execution_mode: self.execution_intent(),
                    backend: SKBackendKind::Faer,
                },
                || state.absorb_batch(&batch.data(), parallelism),
            )?;
            SKStreamDecision::Continue
        }, plan)?;

        // solve at stream end — rank checks ride on G's singular values (≈ σᵢ² of A)
        SKLinearRegressionModel::from_gram(state.solve_with_backend(...)?)
    }
}
}

predict on the model, honest minimum:

#![allow(unused)]
fn main() {
impl<F: SKFloat> SKRegressorPredictor<F> for SKLinearRegressionModel {
    type Error = SKError;
    fn predict<'a, X>(&self, features: X) -> Result<ndarray::Array1<F>, Self::Error>
    where X: TryInto<SKDataView<'a, F>, Error = SKError> {
        let view = features.try_into()?.as_dense()?;
        let (rows, _cols) = view.dim();
        let mut out = ndarray::Array1::zeros(rows);
        azip!((row in view.rows(), value in &mut out) {
            let mut acc = self.intercept_expectation();
            for (coefficient, feature) in self.coefficients.iter().zip(row) {
                acc += coefficient * (*x).get();
            }
            *value = acc;
        });
        Ok(out)
    }
}
}

5. Tests — the numerics-labeled set

#![allow(unused)]
fn main() {
// linear_regression_tests.rs
#[test]
fn recovers_exact_coefficients_to_epsilon() {
    // let x = array![[1.0, 0.0], [0.0, 1.0], [1.0, 1.0], [2.0, 1.0]];
    let true_coefficients = [3.0, -2.0];
    let y: Array1<f64> = x.map(|row| row.dot(&true_coefficients) + 0.5);
    let model = SKLinearRegression::new().build().unwrap().fit(x.view(), y.view()).unwrap();
    assert!((model.coefficients()[0] - 3.0).abs() < 1e-12);
    assert!((model.coefficients()[1] - (-2.0)).abs() < 1e-12);
    assert!((model.intercept() - 0.5).abs() < 1e-12);
}

#[test]
fn rank_deficiency_is_surfaced_not_hidden() {
    // feature columns linearly dependent — rank < p known ahead
    let x = array![[1.0, 2.0], [2.0, 4.0], [3.0, 6.0]];   // second column = 2 · first
    // SVD path: solution & rank reported; NOT silently mis-answered
    let sol = ...fit(x, y) — asserts `rank() == 1` and the minimum-norm answer is returned
}

#[test]
fn streaming_gram_matches_full_array() { /* the streaming state carries exactly
    gram + rhs + row_count; absorb_batch vs full-array G agree to 1e-13 (f64) */
}
}

(walk-through note to keep in the shipped pages: the third test exists to show how the streaming path’s accuracy gate — same fixtures in memory vs streaming — is how normal equations got kept streaming but with the §2 caveat printed as a test contract.)

6. Acceptance walk-through & sources

Chapter sources:

• Golub & Van Loan — §5.3–5.6 (least squares), §12 (SVD applications)
• Trefethen & Bau — lectures 11–12 (SVD), 18–19 (conditioning)
• Higham (2002) — §20 (the condition-number evidence for the normal-equations trap)
• scikit-learn sklearn.linear_model.LinearRegression — the semantics you are matching, and one honest sanity-check: its least-squares dispatch is the SVD-based LAPACK gelsd family — exactly the pseudoinverse path this chapter’s backend code takes

Next: the family where streaming is not an adaptation but the home — gradient descent, state-as-accumulator, and why in-memory training is just “the batch of everything”.

Adding a Streaming Model — SKSGDClassifier and SKSGDRegressor

Pre-read: How to add an Algorithm (§5 planning, §7 streaming driver) and Adding a Linear Model (§2’s condition-number warning — this chapter is the escape route it promised).

This is the chapter where the library’s architecture stops being a compromise. Gradient descent is streaming-native: it consumes one sample (or one batch) at a time, updates an accumulator, and never asks “give me everything”. So the reference code here does not add streaming to a batch algorithm — it presents the estimator as a state machine whose in-memory fit is simply the degenerate case of one batch carrying everything, exactly the “one formulation, two regimes” principle from the Tutorials index.

1. Story & decisions

flowchart TB
    classDef planNode fill:#e0f2fe,stroke:#0284c7,color:#0c4a6e,stroke-width:2px
    classDef computeNode fill:#f0fdf4,stroke:#16a34a,color:#14532d,stroke-width:2px
    classDef modelNode fill:#ede9fe,stroke:#7c3aed,color:#3b0764,stroke-width:2px

    B["SKSGDRegressorBuilder<br>.maximum_batch_samples(1024)<br>.execution_mode(Automatic)"]:::planNode
    B --> E["SKSGDRegressor<br>(immutable estimator)"]:::planNode
    E -->|fit X y| R{"plan"}:::planNode
    R -->|fits in RAM| C["IN-MEMORY<br>one batch = the entire dataset<br>update(batch, parallelism)"]:::computeNode
    R -->|larger than RAM| St["STREAMING<br>sk_run_streaming_driver<br>PartialFit span per batch"]:::computeNode
    C --> Acc["SKSGDState<br>θ, clock, label table (grow-only)<br>θ ← θ − η·∇L"]:::computeNode
    St --> Acc
    Acc --> Ms["SKSGDRegressorModel<br>θ (exportable)"]:::modelNode
    Ms --> P["ŷ = X·θ"]:::modelNode

(Acc marks the accumulator the two regimes share; the classifier carries SKLabelTable in the same state slot.)

2. Theory from zero — gradient descent, without the magic

Both families are supervised, so let us be explicit about the data: you have rows xᵢ ∈ ℝᵖ, and either

• regressor: continuous yᵢ; the loss is squared error L(θ) = ½ Σᵢ (xᵢ·θ − yᵢ)²;
• classifier: labels yᵢ ∈ {0..C−1}; the loss is log-loss over a linear map, with softmax turning the class scores into probabilities, and each class carrying its own weight column.

The property that carries this chapter is convexity. Squared error is a quadratic in θ; log-loss with a linear map is convex. Both are convex, so plain gradient descent converges provably (cited: Bottou 2010 rate analysis; Ruder’s survey as the map). No second-order machinery, no adaptive friend — the bare minimum the modern frameworks demand before their own elaborate bells are even considered.

2.1 The single formula to understand

For a mini-batch Xᵦ (a batch.nrows() × p block) and targets yᵦ, with ridge penalty λ weighed on the coefficient vector:

gradient(θ)   = (1/|B|) · X_Bᵀ (X_B θ − y_B) + λ θ     ← regressor; squared error's gradient
θ_new         = θ − ηₜ · gradient(θ)

Separate facts carry each symbol, and each one is a classic textbook door:

• ηₜ, the learning rate — how far to travel per step. Too high = divergence; too low = days to converge. The two shipped schedules:
  • constant: ηₜ = η₀ — the simplest; tuning sensitivity then lands on the batch length;
  • inverse_scaling: ηₜ = η₀ · t^(−power_t), power_t = 0.25 matching scikit-learn’s default — the Robbins–Monro annealing family, “give a different step scale across the training wave”.
  • optimal (why we analyzed it, not shipped it): scikit-learn’s optimal schedule ηₜ = 1/(α·(t+t₀)) is, in the Pegasos framework, a theorem-coupled rate whose t₀ is derived per Bottou’s practical rules (Bottou 2012, “Stochastic Gradient Descent Tricks”): the constant is dataset-state dependent, and a partial-fit resumable schedule must reproduce that derivation exactly across process restarts to stay a theorem. That is a textbook-level claim worth the chapter’s citation, but not its default: shipping it as pseudo-code without the derivation would hand the reader an unverifiable constant.
• λ (ridge), as the L2 penalty you already know from the PRD — the coefficient shrink is temporal-independent: the same penalty every step; the library holds λ = 0 as legal but noted — a fully quadratic, non-regularized problem is exactly SKLinearRegression, shipped elsewhere by design (convex but closed-form).

Sources:

• Bottou, “Large-scale machine learning with SGD” — the classic: convergence as rates, and the online-vs-batch discussion worth reading before writing
• Ruder, “An overview of gradient descent optimization algorithms”
• Robbins & Monro, “A stochastic approximation method” (1951) — the theorem behind ηₜ ≈ 1/t
• Bottou, “SGD Tricks” (2012) — the practically-correct lemma behind scikit-learn’s optimal
• Shalev-Shwartz et al., “Pegasos: primal estimated sub-gradient solver for SVM” — the theoretical framework tying loss + regularization to convergence with the theoretical learning rate

2.2 Why the classifier’s labels survive streaming (SKLabelTable)

A streaming classifier sees classes as they arrive: batch one shows ["frog", "cat"], batch four introduces "dog" — the label table must grow without changing the meaning of any earlier index. The label table’s deterministic first-occurrence canonicalization is exactly that mechanism: existing labels’ indices are frozen at construction, so batch 4 extends the table from ["frog","cat"] to ["frog","cat","dog"] and index 0/1 are cannot change after first occurrence). A classifier whose coefficient matrix re-indexes labels per batch is buggy in a way that type safety will not catch (shapes stay p × C); §5’s test is the boundary.

3. The combinable state, spelled out

#![allow(unused)]
fn main() {
// crates/sciencekit_linear_model/src/sgd/core_implementation.rs
pub struct SKSGDState {
    /// the schedule's t: number of `absorb_batch` calls the state has survived
    pub clock: u64,
    /// (p, C) — C = 1 for the regressor; f64 state unconditionally, per the anchor
    pub coefficients: ndarray::Array2<f64>,
    pub intercept: ndarray::Array1<f64>,
    /// classifier only — grow-only first-occurrence canonicalization
    pub labels: Option<SKLabelTable>,
}

impl SKSGDState {
    /// THE "partial fit". One batch in; one (θ, t) update out. This is EXACTLY the
    /// closure body handed to `sk_run_streaming_driver` — no outer loop re-derives it.
    pub fn absorb_batch<F: SKFloat>(
        &mut self,
        batch: &ndarray::ArrayView2<F>,
        targets: &ndarray::ArrayView1<'_, f64>,
        parallelism: usize,
        hyper: &SKSGDHyperparameters,
    ) {
        // gradient runs through the vectorized kernels (partial per worker — as in §6
        // of the transformer chapter), so no thread sees another's accumulator:
        let gradient = compute_batch_gradient(batch, targets, &self.coefficients, hyper);
        let rate = self.schedule_rate();   // constant | inverse_scaling by the clock
        self.coefficients = &self.coefficients - &gradient.mapv(|g| rate * g);
        self.intercept = ...;  // same update for the intercept row, per schedule
        self.clock += 1;
    }
}
}

4. Fit, predict, and the twins

SKSupervisedFit goes the same plan-setting four lines as every algorithm (§5 of the anchor) and then, per plan.mode:

#![allow(unused)]
fn main() {
impl<F: SKFloat> SKSupervisedFit<F> for SKSGDRegressor {
    type Model = SKSGDRegressorModel;
    type Error = SKError;
    fn fit<'a, X, T>(&self, features: X, targets: T) -> Result<Self::Model, SKError>
    where
        X: TryInto<SKDataView<'a, F>, Error = SKError>,
        T: TryInto<SKTargetView<'a>, Error = SKError>,
    {
        let view = features.try_into()?.as_dense()?;
        let target = targets.try_into()?.as_continuous()?.into_owned();
        sk_run_operation(
            SKOperationAttributes {
                operation: SKOperationKind::Fit,
                rows: view.nrows(), columns: view.ncols(),
                execution_mode: self.execution_intent(), backend: SKBackendKind::Faer,
            },
            || {
                let mut context = SKExecutionContext::real();
                context.dataset_size_bytes = view.len() as u64 * size_of::<F>() as u64;
                context.dataset_elements = Some(view.len() as u64);
                context.scalar_size_bytes = size_of::<F>() as u64;   // buffer guard precision
                context.access_pattern = SKAccessPattern::Sequential;
                context.grain = SK_ELEMENTWISE_TRANSFORM_GRAIN;   // the update's dominant kernel
                context.batch_size_hint = Some(self.maximum_batch_samples());
                let plan = sk_resolve_execution_plan(self.execution_intent(), &context)?;

                use SKExecutionMode::*;
                match plan.mode {
                    InProcessSynchronous | InProcessAsynchronous => {
                        // THE degenerate batch: the entire array is one batch
                        let mut state = SKSGDState::zeros(view.ncols(), self.hyperparameters());
                        state.absorb_batch(&view, &target.view(), plan.parallelism, &self.hyperparameters());
                        Ok(SKSGDRegressorModel::from_state(state))
                    }
                    OutOfCoreMemoryMapped => Err(SKError::ExecutionModeIncompatible {
                        mode: "OutOfCoreMemoryMapped", pattern: "sequential-batch SGD",
                    }),
                    OutOfCoreStreaming => Err(SKError::ExecutionModeIncompatible {
                        mode: "OutOfCoreStreaming", pattern: "in-memory fit entry point (use fit_streaming_source)",
                    }),
                    Automatic => unreachable!("resolve always returns a concrete mode"),
                }
            },
        )
    }
}
}

and the streaming twin (the natural home):

#![allow(unused)]
fn main() {
impl SKSGDRegressor {
    pub fn fit_streaming_source<S, F>(
        &self, source: &mut S, plan: &SKExecutionPlan,
    ) -> Result<SKSGDRegressorModel, SKError>
    where
        F: SKFloat + Send,
        S: SKLazySource<F> + Send,
    {
        let mut state = SKSGDState::zeros(...);
        sk_run_streaming_driver(source, &mut state, |batch, parallelism, state| {
            sk_run_operation(
                SKOperationAttributes {
                    operation: SKOperationKind::PartialFit,   // one span per batch
                    rows: batch.data().nrows(), columns: batch.data().ncols(),
                    execution_mode: self.execution_intent(), backend: SKBackendKind::Faer,
                },
                || {
                    state.absorb_batch(&batch.data(), &batch_targets, parallelism, &self.hyperparameters());
                    SKStreamDecision::Continue   // keep streaming until the source ends
                },
            )?,
        }, plan)?;
        Ok(SKSGDRegressorModel::from_state(state))
    }
}
}

The exact update closure (sketch carried faithfully into chapter 3):

one batch in: Xᵦ, yᵦ → gradient(θ) ← (1/|B|) (Xᵦ θ − yᵦ)ᵀ Xᵦ + λ θ
              rate ← schedule(η₀, t)         # one of {constant, inverse_scaling}
              θ   ← θ − rate·gradient        # the ONLY sequential-mutating side of the family
              t   ← t + 1

5.6 predict and the label table’s reverse map

The classifier’s predict reduces probabilities across categories via the same label table — argmax over the class weights, then SKLabelTable::label_of(index) decodes the integer back into its frozen label string. The regressor simply returns one F per row — matching the anchor §3.2 rule: SKRegressorPredictor<F> returns Array1<F>.

The structural rule this encodes: label lookups always go through the table, never an ad-hoc map in the model — §5’s grow-only test is the enforcement.

5. Tests — the streaming-native four

#![allow(unused)]
fn main() {
// sgd_tests.rs
#[test]
fn in_memory_equals_single_batch_streaming() {
    // identical answers: fit(whole array) == fit_streaming(one batch holding the array)
    // — the anchor principle, tested at the numerical level
}

#[test]
fn learning_rate_schedule_matches_specified() {
    // schedule-only change: constant vs inverse_scaling on the SAME fixture batch
    // (θ paths differ *by the schedule formula*, nothing else)
}

#[test]
fn labels_table_is_grow_only_in_stream() {
    // batch 1: ["frog","cat"]; batch 4 introduces "dog"
    // asserts: pre-"dog" indices unchanged after the fourth batch, first-occurrence
    // — inference at `sk_canonicalize_labels` — and prediction decodes via `label_of`
}

#[test]
fn per_batch_partial_fit_spans() {
    // asserts count of PartialFit spans == count of batches — the legible trace
    // that makes streaming jobs debuggable rather than one opaque fit span
}
}

6. Acceptance & sources — mapped, not repeated

Standing study sources (each also explained §2):

• Bottou, “Large-scale machine learning with SGD” (2010) & “SGD Tricks” (2012)
• Ruder (see §2)
• Shalev-Shwartz et al., Pegasos (see §2)
• Robbins & Monro (see §2)
• scikit-learn sklearn.linear_model.SGDClassifier/SGDRegressor documentation and the partial_fit protocol docs — the semantics referenced in this chapter, including the honest caveat on optimal.

Next: the family where streaming does not apply — and explaining why is the lesson.

Adding a Neighbors Model — SKKNeighborsClassifier and SKKNeighborsRegressor

Pre-read: How to add an Algorithm (§5 access patterns, §8 SKMappableSource)

This is the chapter that keeps the book honest. Every algorithm so far followed one shape: fit computes, predict applies. k-NN breaks the shape — and understanding why it cannot be repaired is the lesson of the chapter. In sciencekit terms:

• fit does almost nothing: it stores a view of the training base (or the path/SKMappableSource handle for the memory-mapped case).
• predict does everything: it needs, for a query row, the k nearest rows of the base — which is random access to any stored row at any time.

A sequential streaming path (OutOfCoreStreaming) therefore does not fit this family, no matter how it is dressed up. The formulation that fits k-NN’s pattern is SKAccessPattern::RandomAccess with an SKMappableSource base: the out-of-core version is the memory-mapped base, not a streaming combiner.

1. Story & decisions

flowchart TB
    classDef fitNode fill:#e0f2fe,stroke:#0284c7,color:#0c4a6e,stroke-width:2px
    classDef predictNode fill:#f0fdf4,stroke:#16a34a,color:#14532d,stroke-width:2px
    classDef coreNode fill:#ede9fe,stroke:#7c3aed,color:#3b0764,stroke-width:2px

    F["SKKNeighborsClassifierBuilder<br>.nearest_neighbors_count(5)<br>.weights(Uniform | Distance)<br>.execution_mode(Automatic)"]:::fitNode
    F --> Fit["fit(X, y)<br>· stores base rows as SKMappableSource<br>(in RAM or memory-mapped)<br>· freezes SKLabelTable (grow-only)"]:::fitNode
    Fit --> Q["predict(X_query)"]:::predictNode
    Q --> D["for each query row:<br>sk_squared_euclidean_distance_matrix<br>chunked over base rows (cache grain)"]:::predictNode
    D --> T["top-k via argpartition over per-base chunks<br>— O(n) expected, not O(n log n) sort"]:::predictNode
    T --> W["uniform vote | inverse-distance weights"]:::predictNode
    W --> O["majority via SKLabelTable<br>(classifier) | mean of neighbors (regressor)"]:::coreNode

(the diagram is fit heavy: the predict spine is where all compute lives — the shape that distinguishes this family).

2. Theory from zero — neighbor search

The problem. Given a base of named rows (with labels or target values) and a query row x, return the k stored rows with the smallest distance to x, then reduce those k into one answer (majority label for a classifier; mean for a regressor).

To build the distance: for base row b and query row x, ‖x − b‖₂² = x·x + b·b − 2·x·b, with two squared norms precomputable per row — so the cost of a row pair reduces to the dot product x·b, which is exactly what a GEMM-shaped operation summation can take down to the library’s grain. That is to say: the whole distance structuring of brute-force k-NN is a matrix-multiplication-shaped question, and this is the reason the library’s pairwise kernels (sk_squared_euclidean_distance_matrix, sk_euclidean_distance_matrix, sk_manhattan_distance_matrix, sk_cosine_distance_matrix) live in sciencekit_math::pairwise with that shape already.

2.1 The puzzle of finding “the nearest few”

Naively: compute all base distances for the query, sort them, take the first k — O(n log n). Selection theory rides in and pays rent immediately: you do not need a sorted order, you need the k smallest — the answer to “subset selection under an order” is argpartition-structured quickselect at O(n) expected cost (with a small constant in practice). And at a base of 10⁷ rows and 10⁵ queries, the “sort everything once upfront” precomputation path buys the code nothing but memory pressure, because the query side is unknown at fit time.

2.2 The curse of dimensionality — your honest adversary

Every neighbor-model chapter carries two classic facts and one modern nuance. Classic:

• The curse: in a d-dimensional unit cube, the volume the shell of width ε concentrates almost all the mass for large d; equivalently, the distance between two random points becomes concentrated — so “the nearest” and “any random point” go asymptotically indistinguishable as d grows.
• The kd-tree’s territory: binary space partitioning makes queries O(log n) when log-dim branch-cutting actually prunes candidates, which degrades to near-brute-force in high dimensions. The practical crossover sits low: past roughly ~10–15 dimensions, trees lose to brute force on most realistic data.

The modern nuance this chapter also hands you (the one that makes the decision table non-trivial): metrics this library knows — cosine, manhattan — are metric-dependent in their search-structure costs, and the index raises and falls on that. The distance is not a “toggle” in the estimator; it is a design choice with structure-level consequences. This is exactly the kind of non-obvious decision this chapter exists to flag, and carries citations to read the source before choosing:

Sources:

• Bentley, “Multidimensional binary search trees used for associative searching” (1975) — the kd-tree paper
• Cover & Hart, “Nearest neighbor pattern classification” (1967) — the 1-NN ≤ 2×Bayes bound; note the k→∞, k/n→0 consistency result is Stone (1977) — often misattributed, cite Stone
• Friedman, “An improved algorithm finding nearest neighbors” — improved kd-tree search ordering; the separate best-bin-first variant is Beis & Lowe (1997) — cite both
• Muja & Lowe, “Fast approximate nearest neighbors with automatic algorithm configuration” (FLANN) — the practical decision territory of modern neighbor search
• Wikipedia: k-d tree and Nearest neighbor search — the map of exact vs approximate avenues
• scikit-learn sklearn.neighbors — KNeighborsClassifier/KNeighborsRegressor semantics (weights, algorithm ‘auto’ — brute/kd/ball — and the exact behavior of distance weights)

3. The model’s shape

#![allow(unused)]
fn main() {
// crates/sciencekit_neighbors/src/k_neighbors/core_implementation.rs
pub enum SKNeighborWeights { Uniform, Distance }

pub struct SKKNeighborsModel<F: SKFloat, S: SKMappableSource<F>> {
    base: S,                        // the training rows — the whole storage question
    targets: ndarray::Array1<F>,  // owned continuous targets (regressor only)
    labels: Option<SKLabelTable>,   // classifier only — frozen, grow-only
    nearest_neighbors_count: usize,
    weights: SKNeighborWeights,
    /// the intent recorded from the builder at fit time — the model's invoke of
    /// execution_intent() reads it, while *resolution* still happens per operation
    /// against a fresh SKExecutionContext (anchor §4.1)
    execution_intent: SKExecutionMode,
}

impl<F: SKFloat, S: SKMappableSource<F>> SKKNeighborsModel<F, S> {
    pub fn execution_intent(&self) -> SKExecutionMode { self.execution_intent }
}
}

Two notes before the code:

1. The model owns an SKMappableSource for the base. In memory this is a borrowed SKDataView asserting row-major contiguous layout (via the zero-copy conversions chapter’s rule: the view is not copied); for the memory-mapped regime, a memmap-backed implementation hands row(index) out from the file — same trait, same algorithm code.
2. predict runs per query row, and the base row iteration is chunked over the cache grain — which is exactly the grain discipline from the anchor chapter (§6): the chunk size, not the rows’ total, is what carries the parallelism.

4. The crown code

#![allow(unused)]
fn main() {
// Regressor: responses preserve the model scalar. The classifier sibling
// implements SKClassifierPredictor<F> instead, returning Array1<i64> labels.
impl<F, S> SKRegressorPredictor<F> for SKKNeighborsModel<F, S>
where
    F: SKFloat,
    S: SKMappableSource<F> + Sync,
{
    type Error = SKError;
    fn predict<'a, X>(&self, queries: X) -> Result<ndarray::Array1<F>, Self::Error>
    where X: TryInto<SKDataView<'a, F>, Error = SKError> {
        let queries = queries.try_into()?.as_dense()?;
        sk_run_operation(
            SKOperationAttributes {
                operation: SKOperationKind::Predict,
                rows: queries.nrows(), columns: queries.ncols(),
                execution_mode: self.execution_intent(), backend: SKBackendKind::Faer,
            },
            || {
                let mut context = SKExecutionContext::real();
                // the plan rides from the estimator's SKBuilderState (the model owns no
                // plan — anchor §4.1); resolution happens per predict
                context.dataset_size_bytes = queries.len() as u64 * size_of::<F>() as u64;
                context.dataset_elements = Some(queries.len() as u64);
                context.access_pattern = SKAccessPattern::RandomAccess;   // the family's pattern
                context.grain = SK_BINARY_COMBINE_GRAIN; // pairwise/candidate kernels;
                // (a dedicated GEMM grain is not yet measured — revisit once one is measured)
                let plan = sk_resolve_execution_plan(self.execution_intent(), &context)?;

                let rows = queries.nrows();
                let mut out = ndarray::Array1::<F>::zeros(rows);

                // A per-query answer: every query row depends only on the base,
                // so there is no reduction and no shared accumulator — par_azip! is safe.
                if plan.parallelism > 1 {
                    par_azip!((query in queries.rows(), out in &mut out)) {
                        *out = self.answer_for_query(&query, plan.parallelism);
                    }
                } else {
                    azip!((query in queries.rows(), out in &mut out)) {
                        *out = self.answer_for_query(&query, 1);
                    }
                }
                Ok(out)
            },
        )
    }
}

// One query row: distance to every base row, then top-k selection, then reduce.
impl<F: SKFloat, S: SKMappableSource<F>> SKKNeighborsModel<F, S> {
    fn answer_for_query(
        &self,
        query: &ndarray::ArrayView1<F>,
        parallelism: usize,
    ) -> F {
        // squared distances: x·x + b·b − 2·x·b, per-chunk over the mappable base
        let distances = sciencekit_math::pairwise::sk_squared_euclidean_distance_matrix(
            query.insert_axis(ndarray::Axis(0)),   // 1 × p
            self.base_matrix_view(),               // base rows, chunked from the source
            parallelism,
        );
        // top-k: partition-select the k smallest — O(n) expected, not O(n log n) sort
        let candidates = top_k_selection(&distances.row(0), self.nearest_neighbors_count);
        self.answer_for(&candidates, &self.weights)
    }
}
}

(prose: top_k_selection-shaped partitioning carries through ndarray partitioning; the test in §5 walks the exact selection kernel naming — argpartition shape as the library’s helper name once it lands in sciencekit_math::pairwise — and the crown code uses the same SKStreamDecision-free, pure-kernel shape as every other regime; the chapter’s code here is the walk-through of the dual-form principle applied to the random-access family.)

The two answer rules:

# uniform:  classifier — the mode of the k neighbors' labels (ties → the label with
#           the smallest canonical index; deterministic even when tied)
#           regressor  — the mean of the k neighbors' targets

# distance: both families weight each neighbor by 1 / dist, then normalize. The reviewed
#           edge is the zero-distance neighbor (dist = 0): scikit-learn assigns it weight
#           1.0 and then normalizes — that is the rule implemented here; the ε-clamp is a
#           *separate* documented safety choice, not "the same decision" as scikit-learn

5. Tests — what can honestly be asserted here

#![allow(unused)]
fn main() {
// k_neighbors_tests.rs
#[test]
fn count_of_neighbors_respects_base_size_at_construction_time() {
    // builder validation: nearest_neighbors_count ≤ base rows is evaluated at fit,
    // because the base's shape is data, not configuration — the error taxonomy's
    // InvalidHyperparameter is surfaced AT fit instead of double-checking downstream
}

#[test]
fn classifier_majority_ties_deterministic() {
    // k=4, base labels [cat, dog, cat, dog] — tie → first-occurrence label wins:
    // SKLabelTable's determinism is THE answer, not a coin toss
}

#[test]
fn regressor_row_lies_along_neighbors() {
    // squared-distance reduction: best-k rows hold the winner; "distance" weights
    // and the ε guard are asserted with a duplicate base row (dist = 0) — the §NN
    // guard question shipped as tests.
}

#[test]
fn memory_mapped_base_matches_in_memory_base() {
    // same base data, one via RAM view, one via memmap-backed SKMappableSource —
    // predictions are exact-equal; the third circuit rides through the trait seam
}
}

6. Acceptance & sources

Chapter sources (see §2.1’s citations): Bentley, Cover & Hart, Muja & Lowe, the sklearn.neighbors docs and its npy-neighbour algorithm notes on the algorithm = "auto" contract.

7. What the four chapters added up to

You have now written four families: two state-based transformers with streaming, a dense linear model with solver strategy, an SGD-family with an accumulator, and this chapter’s random-access exception. Look back: every one of them had identical skeleton files (builder, estimator, model, traits, plan wiring, observability, tests) and different numeric cores. That articulation — shell the same, math the different — is sciencekit’s thesis, and it is your job as a contributor to hold the shell invariant while freely inventing the math inside.

Contribute

Where help is needed

The roadmap is large by design — every algorithm is an isolated, well-scoped opportunity:

Open an issue or grab one labeled good first issue — specs are written before code, so you always start with clear requirements.

How we work

1. Spec first (OpenSpec): each change gets a proposal, spec deltas and tasks on a dedicated -openspec branch.
2. TDD: tests are written first and must fail; implementation follows minimally.
3. Small PRs: every change alters the minimum necessary; work happens in git worktrees, never directly on main.
4. Independent review: an independent reviewer validates the spec was met before merge.
5. Acceptance: every implementation runs with lots and little data, under concurrency, exports its model and produces metrics.

Full engineering rules live in AGENTS.md and the product truth in the PRD.

Getting started

git clone https://github.com/neurono-ml/sciencekit.git
cd sciencekit
cargo build && cargo test    # Rust 1.85 pinned via rust-toolchain.toml

License

Apache-2.0 — including patent protection. By contributing you agree your contributions are licensed under it.

Documentation Guide

Instructions for agents and contributors writing or updating this book. Follow them so every page keeps the same visual language and quality bar.

Where things live

Branch workflow (mandatory)

• The documentation lives on branch docs/documentations — never commit directly to it from main.
• Every docs change happens in a worktree: temporary/worktrees/docs/<change-name> on branch docs/<change-name>, merged into docs/documentations when done.
• Pushing to docs/documentations automatically rebuilds and deploys the site.

Local preview

mdbook serve docs --open     # live reload at http://localhost:3000
mdbook build docs            # output in docs/book/ (never committed)

Language rules

• The book is written in English. Everything durable in this repository is English.
• The PRD (docs/PRD.md) is the source of truth; translate its concepts faithfully, never invent scope.
• Code identifiers stay exactly as they are (SKStandardScaler, sk_train_test_split) — do not translate names.

Component catalog

Use the ready-made components instead of raw styling. All accept standard HTML attributes.

Hero + announcement pill

<div class="sk-hero reveal">
  <div class="sk-announce">
    <span class="sk-announce__tag">New</span>
    Supporting text
  </div>
  <h1 class="sk-hero__title">Big statement title</h1>
  <p class="sk-hero__subtitle">Supporting sentence.</p>
</div>

Cards grid

<div class="sk-cards reveal">
  <div class="sk-card">
    <span class="sk-card__icon sk-icon--sky">⚡</span>
    <div class="sk-card__title">Title</div>
    <p>Description.</p>
  </div>
</div>

Icon tints: sk-icon--sky | violet | emerald | amber | rose | neutral. Add sk-card--flow for gradient-washed feature cards.

Callout boxes

<div class="sk-box sk-box--info">…</div>

Variants: sk-box--info (blue), sk-box--tip (violet), sk-box--success (green), sk-box--warning (amber), sk-box--danger (red).

Stats band

<div class="sk-stats reveal">
  <div class="sk-stats__item"><div class="sk-stats__value js-count" data-target="39" data-suffix="+">39+</div><span class="sk-stats__label">algorithms mapped</span></div>
</div>

js-count animates the number when it scrolls into view.

Pills (status badges)

<span class="sk-pill sk-pill--emerald">Done</span>

Variants: sk-pill--sky | violet | emerald | amber | rose | neutral.

Progress bar

<div class="sk-progress"><div class="sk-progress__bar sk-progress__bar--amber" style="width:18%"></div></div>
<div class="sk-progress__meta"><span>label</span><span>18%</span></div>

Bar colors mirror pills: --sky, --violet, --emerald, --amber.

Timeline (roadmap)

<div class="sk-timeline">
  <div class="sk-timeline__item is-active reveal">
    <div class="sk-timeline__card">
      <div class="sk-timeline__head">
        <span class="sk-timeline__title">Phase N — Name</span>
        <span class="sk-pill sk-pill--amber">In progress</span>
      </div>
      <p>Summary.</p>
    </div>
  </div>
</div>

Item states: is-done (green dot), is-active (amber dot), default (violet dot).

Donut chart & CSS bars

Copy the SVG structure from Algorithms and recompute segments: circumference = 2π × 80 ≈ 502.65; segment = count / total × 502.65; each circle’s stroke-dashoffset is the negative sum of all previous segments. Keep legend counts in sync.

For horizontal bars reuse .sk-bars rows as seen in Plan & Roadmap.

Terminal window

<div class="sk-window">
  <div class="sk-window__bar">
    <span class="sk-window__dot sk-window__dot--r"></span>
    <span class="sk-window__dot sk-window__dot--y"></span>
    <span class="sk-window__dot sk-window__dot--g"></span>
    <span class="sk-window__title">bash</span>
  </div>

  ```bash
  cargo test --workspace
  ```

</div>

Statement panel & closing CTA

<div class="sk-statement reveal"><p class="sk-statement__text">One strong sentence.</p></div>

<div class="sk-cta reveal">
  <div class="sk-cta__title">Call to action</div>
  <p>Why act.</p>
  <a class="sk-btn sk-btn--star" href="#">⭐ Label</a>
</div>

Buttons

<a class="sk-btn sk-btn--primary" href="./page.html">Primary action</a>       <!-- dark -->
<a class="sk-btn sk-btn--secondary" href="#">Secondary</a>                    <!-- outline -->
<a class="sk-btn sk-btn--star" href="...">⭐ Star-flavored</a>                <!-- amber -->

Content rules

1. Add new chapters to SUMMARY.md under the right part (Project / Engineering / Community / Development Explained).
2. Keep numbers honest: algorithm counts, phase progress and status pills must match the roadmap state.
3. Prefer components over custom inline CSS; extend skin/custom.css only for genuinely reusable pieces.
4. Add the reveal class to major blocks for scroll animation — sparingly.
5. Internal links use relative paths ending in .html (./algorithms.html).
6. Run mdbook build docs and fix every warning before pushing; create-missing = false means broken links fail the build.
7. Diagrams must be Mermaid, never ASCII. Write a ```mermaid fence; the book renders it client-side (see AGENTS.md → Documentation). Add an accTitle/accDescr to every diagram so it is accessible.
8. Mermaid colors must use the fixed --sk-* pairs. Every classDef/style declares fill + stroke + color from the approved pairs below — never fill alone, so text stays legible in both book themes:
   • blue: fill:#e0f2fe,stroke:#0284c7,color:#0c4a6e
   • green: fill:#f0fdf4,stroke:#16a34a,color:#14532d
   • orange: fill:#fff7ed,stroke:#ea580c,color:#7c2d12
   • violet: fill:#ede9fe,stroke:#7c3aed,color:#3b0764
   • amber: fill:#fef3c7,stroke:#d97706,color:#78350f
   • rose: fill:#ffe4e6,stroke:#e11d48,color:#7f1d1d
9. Mermaid labels must not clip. Keep each line at or below ~24 characters, break long node labels with explicit <br/>, and put spaces around separators in edge labels (|InProcess* | whole array|, not |InProcess*|whole array|). Unbreakable tokens (type names, paths) go on their own <br/> line.
10. Equations must be LaTeX, not plain text. Write $...$ (inline) or $$...$$ (display) so KaTeX renders them. Do not paste ‖a‖², Σ, √ etc. as plain text.
11. Every image needs descriptive alt text. Diagrams carry accTitle/accDescr; images (.png/.svg) must have an alt attribute that explains the content as read, not just a filename or raw code.

Checklist for a docs change

• Worktree created on temporary/worktrees/docs/<name>
• Content in English, consistent with the PRD
• New pages registered in SUMMARY.md
• Components used per catalog above
• mdbook build docs passes with zero warnings
• Merged into docs/documentations (Pages deploys automatically)

Developer Book

Welcome to the Developer Book — a textbook inside the sciencekit documentation. The other pages tell you what sciencekit is. These chapters tell you how sciencekit was built — decision by decision, pull request by pull request — and, more importantly, why each choice was made, and what alternatives were deliberately thrown away.

You do not need a background in Rust or in numerical computing to read it. You need:

• Basic programming — you know what a function, a loop, a variable and an if are.
• Curiosity — you want to know why, not just what.
• Patience — some ideas (traits, generics, zero-copy) take a chapter or two to click.

How to read this book

Each chapter corresponds to one real pull request that changed the code, and each is written as a self-contained lesson with the same skeleton:

1. The problem we were solving — a plain-language story, before any math.
2. The decisions — and the roads not taken — the choices, and the alternatives we rejected.
3. The concepts, taught from zero — the Rust and numerical ideas the chapter needs.
4. Each object, explained — every public struct, trait, enum and function, and why it exists.
5. A real walk-through, using the tests — the actual tests from the repository, read line by line.
6. Look inside — a diagram and a strong insight.
7. Recap — what you should take away.

Why we use the tests as examples

sciencekit is built with test-driven development: the tests in the *_tests.rs companion modules are the executable specification of every function. When a chapter says “here is how this works”, it points you at the real test that proves it. You can read the code, then run cargo test and watch it pass. The test is the proof; the chapter is the explanation.

The order matters

The chapters follow the order the code was written. Earlier chapters build concepts that later chapters assume. In particular, Core contracts is the keystone — almost every later chapter builds on the types it introduces. If you only read one chapter first, read that one.

Bootstrap: workspace & CI     — the skeleton
Core contracts                — the vocabulary (traits, errors, views)  ← keystone
Math kernels & BLAS           — the numeric substrate
Execution & observability
The host-centric backend
Module hygiene                — how the code is organised as it grows
Parallel kernels & streaming

How to keep this book honest

sciencekit is still being built. When new algorithms land, or a bug fix changes behaviour, these chapters must change too. The final chapter, Maintaining this book, explains that rule: new features get a new chapter; fixes and refactors update an existing one. If you find a chapter that no longer matches the code, that is a bug — in the book, not in the code — and it should be fixed the same way a code bug is.

Bootstrap: the Cargo workspace & CI

Every machine-learning library is, underneath, a collection of many small pieces that must agree on shared rules: how numbers are typed, how errors are reported, how code is formatted, and how it is proven to work. When sciencekit started, it was an empty repository — no code at all. The very first engineering decision was therefore not which algorithm to write, but how to lay the ground so that every later chapter has a stable place to stand. This chapter is about that foundation: a Cargo workspace that splits the library into sub-crates, a pinned Rust toolchain, and an automated CI pipeline that guards the whole thing.

1. The problem we were solving

Imagine you are going to build a house. You do not start by picking paint colors. You start with the foundation, the plumbing, the frame — the parts that are invisible once the house is finished but that everything else depends on.

A Rust library has the same hidden first step. Before you can write a single estimator, you need to answer boring-but-critical questions:

• Where does the code live? One giant folder, or several small, related folders?
• Which version of the Rust language are we writing for? Rust evolves every six weeks; we must say exactly which one we target.
• How do we prove the code works? Running the tests by hand on your laptop is not enough. A machine must re-run them automatically on every change, forever, so a mistake is caught the moment it is made.

If you skip these questions, you get a pile of code that:

• compiles on the author’s machine but not anyone else’s;
• mixes unrelated concerns in one huge, slow-to-compile crate;
• silently breaks because nobody re-ran the tests.

In sciencekit’s case the very first commit had no code at all — the repository was freshly initialized under the name neurono-ml. So the first task, called Phase 0.1 “bootstrap”, was to create the skeleton: a workspace, a pinned toolchain, and a CI pipeline. Nothing else. Two pull requests did this:

• PR 96b77ab — “bootstrap sciencekit workspace”: added Cargo.toml (the workspace root), rust-toolchain.toml (the pinned toolchain), a LICENSE, and a nearly-empty src/lib.rs.
• PR dd3b524 — “add CI workflow”: added .github/workflows/ci.yml with five automated gates.

Everything you will read in later chapters — sciencekit_common, sciencekit_math, every algorithm — sits on top of this skeleton. That is why the chapter is called “bootstrap”: a library cannot build itself up until this scaffolding exists.

2. The decisions — and the roads not taken

A bootstrap is almost entirely decisions about structure. Here are the four big ones, and the alternatives we consciously rejected.

Two smaller but telling decisions in the same spirit:

• License: Apache-2.0 — chosen for the project (visible in Cargo.toml), standard for an international open-source library.
• src/lib.rs declares #![forbid(unsafe_code)] — a library this concerned with correctness opts out of unsafe Rust by default, so no unsound code can sneak in silently. (A later phase reintroduces a controlled, reviewable amount of unsafe for SIMD — but only behind explicit, approved features.)

3. The concepts, taught from zero

Before we read the files, we need three ideas from scratch. If you already know them, skim; if not, this is the whole foundation.

3.1 What is a “crate”? What is a “workspace”?

A crate is Rust’s unit of compilation and distribution. Think of it as a package: it has a name, a version, a set of dependencies, and a set of source files. When you build, Rust compiles each crate as a unit.

A workspace is a collection of crates that are developed together in one repository. They share a single Cargo.toml at the root that lists them. You can build, test, and release them all with one command (cargo build --workspace).

Why does this matter for a library? Three reasons:

1. Isolation. If sciencekit_math (linear algebra) pulls in a heavy dependency, that dependency only affects sciencekit_math. The other crates do not see it, do not recompile because of it, and are not slowed down by it.
2. Small PRs. Each algorithm area becomes its own crate, added in its own change. Reviewers read one focused change instead of a megamorphic one.
3. Reuse. Sub-crates can depend on each other cleanly — sciencekit_math depends on sciencekit_common, for instance — using path dependencies (relative paths inside the workspace) rather than downloading from the internet.

Think of a workspace like a monorepo: one repository holding several related projects that know about each other.

3.2 Edition 2024 and MSRV

Rust has two “versions” knobs that beginners often confuse:

• The edition is a language compatibility mode. Every few years the Rust team makes a release that may change how some syntax behaves. To avoid breaking all existing code, they ship these as “editions” (2015, 2018, 2021, and now 2024). Picking edition = "2024" just says “write this crate using the 2024 rules of the language.” It is set once, in Cargo.toml.

• MSRV stands for Minimum Supported Rust Version. It is the oldest Rust compiler that this project promises to work with. Our MSRV is Rust 1.85 (set via rust-version = "1.85" in Cargo.toml). This is a compatibility promise to users: they can build with 1.85 or anything newer.

There is also the subtle difference between the two places we record the toolchain:

• rust-version = "1.85" in Cargo.toml is a promise to the world — it says “1.85 is the minimum you need.” It is metadata.
• rust-toolchain.toml is a declaration for our own machines — it says “use exactly 1.85, no matter what is installed.” It enforces reproducibility for contributors and CI.

The difference matters: one is a guarantee, the other is a lock.

3.3 What is a CI pipeline? What are “gates”?

CI stands for Continuous Integration. The idea: instead of hoping everyone runs the tests locally, a remote server automatically runs a checklist of checks on every code change. If any check fails, the change is blocked — it cannot be merged.

Each check is a gate: a rule the code must pass to be accepted. Think of it like airport security — you do not board until every checkpoint clears. The sciencekit CI has five gates, and each guards a different failure:

The msrv gate is the subtle one: our own machines might run Rust 1.90, but if someone on 1.85 is promised support, the code must actually build there. A dedicated job running exactly 1.85 proves it. That is the “reproducibility promise” made real.

3.4 The naming convention: SK / sk_

One more foundation stone, defined in PRD §3.4, will be visible in every later chapter, so it is worth learning now:

• No abbreviations, except a single project prefix. Names must be complete and self-explanatory: maximum_number_of_iterations, not max_iter; nearest_neighbors_count, not k.
• The project prefix sk / SK is the only allowed abbreviation.
• Structs and traits get SK + PascalCase: SKEstimator, SKStandardScaler, SKExecutionMode.
• Free functions, variables, and modules get sk_ + snake_case: sk_train_test_split, sk_resolve_execution_plan.
• Methods (functions inside an impl block) get no prefix: fit, predict, execution_mode.
• Crates always keep the full name sciencekit_* — never abbreviated.

Why? Machine-learning code is full of cryptic one-letter names (k, c, X, y). The project deliberately bans them for readability and self-documentation. The sk prefix is a namespace marker: it tells readers “this object belongs to sciencekit” the way a brand marks its products.

4. Each object, explained

For the bootstrap, the “objects” are configuration files and the workflow, not Rust structs. Each one is explained below in the same mini-structure.

4.1 Cargo.toml (workspace root)

What it is. The workspace’s identity card and command center. It declares the crate sciencekit and the workspace that will hold the sub-crates.

Why it exists. It is where the workspace is born. It records the name, version, edition, license, MSRV, and — crucially — the empty list of workspace members. Here is the real content (from PR 96b77ab):

#![allow(unused)]
fn main() {
[package]
name = "sciencekit"
version = "0.1.0"
edition = "2024"
license = "Apache-2.0"
rust-version = "1.85"
description = "Machine Learning library in Rust reimplementing all of scikit-learn — extreme performance, zero-copy and native out-of-core."
repository = "https://github.com/neurono-ml/sciencekit"

[workspace]
resolver = "3"
members = []
}

What we rejected. Filling members with a dozen not-yet-written sub-crates. The decision (PRD §3.1) is to add sub-crates one per change, so members starts empty and grows deliberately. The resolver = "3" line picks the modern dependency-resolution strategy (required for edition 2024 workspaces).

4.2 rust-toolchain.toml

What it is. The reproducibility lock. It tells any machine that opens this repository, “install and use exactly this Rust toolchain.”

Why it exists. To guarantee that a clean clone compiles identically for every contributor and for CI. The real content:

# Pin the exact Rust toolchain for reproducible builds (PRD §3.6).
# A clean clone compiles with exactly 1.85, edition 2024, regardless of the
# contributor's installed Rust. Upgrading the channel is an explicit, reviewable change.

[toolchain]
channel = "1.85"
components = ["rustfmt", "clippy"]
profile = "minimal"

Note the components line: it asks for rustfmt (the formatter) and clippy (the linter) to be installed alongside the compiler — the two tools our CI gates depend on. profile = "minimal" avoids installing extra docs/components we do not need, keeping clones fast.

What we rejected. Not pinning, or pinning to a moving “stable” channel. Both would break reproducibility.

4.3 src/lib.rs (the umbrella crate)

What it is. The root source file of the sciencekit crate — currently nearly empty.

Why it exists. It is the future umbrella: PRD §3.1 and §3.3 describe sciencekit as the crate that re-exports the public surface of all the sciencekit_* sub-crates, so users can write use sciencekit::prelude::* and get everything. At bootstrap there is nothing to re-export yet, so it carries only the two attribute declarations:

#![allow(unused)]
fn main() {
//! `sciencekit` — Machine Learning library in Rust reimplementing all of
//! scikit-learn, with extreme performance, zero-copy and native out-of-core.
//!
//! This crate is the umbrella re-export of the `sciencekit_*` sub-crates.

#![forbid(unsafe_code)]
#![warn(missing_docs)]
}

#![forbid(unsafe_code)] refuses any unsafe block (a correctness stance). #![warn(missing_docs)] warns when public items lack documentation (a quality gate — the whole library is to be documented).

What we rejected. Putting the umbrella crate logic inside the workspace root crate from day one. It will be a separate crate under crates/ once there is something to re-export.

4.4 .github/workflows/ci.yml (the CI pipeline)

What it is. The automation recipe. It is a GitHub Actions workflow: a YAML file that tells GitHub’s servers what to run on every push and every pull request.

Why it exists. To make the five gates (Section 3.3) run automatically and block bad changes. It triggers on push to main and on every pull_request. Let us look at the shape of two jobs, then the test job, then the MSRV job.

The fmt job checks formatting:

  fmt:
    name: Formatting (cargo fmt)
    runs-on: ubuntu-latest
    steps:
      - name: Install Rust 1.85 toolchain (rustfmt)
        uses: dtolnay/rust-toolchain@1.85
        with:
          components: rustfmt
      - name: cargo fmt --all -- --check
        run: cargo fmt --all -- --check

Note the trick in the last step: cargo fmt --check does not reformat the code; it only reports whether the code already matches the canonical format. If a contributor forgot to format, this job fails.

The clippy job runs the linter with warnings promoted to hard errors:

      - name: cargo clippy --workspace --all-targets -- -D warnings
        run: cargo clippy --workspace --all-targets -- -D warnings

The -D warnings flag is the teeth: a mere warning (e.g. an unnecessary .clone()) becomes a build error. This is deliberate strictness — the project would rather block a merge than ship a smell.

The test job runs everything:

      - name: cargo test --workspace --all-targets
        run: cargo test --workspace --all-targets

--workspace builds and tests every sub-crate; --all-targets also builds test code, examples, and benchmarks so that their compilation is checked too.

And the msrv job re-runs build + test on the exact promised minimum, closing the reproducibility loop:

  msrv:
    name: MSRV gate (exact 1.85)
    runs-on: ubuntu-latest
    steps:
      - name: Install exact Rust 1.85 toolchain (rustfmt + clippy)
        uses: dtolnay/rust-toolchain@1.85
        with:
          components: rustfmt, clippy
      - name: cargo build + test on exact 1.85
        run: |
          cargo build --workspace --all-targets
          cargo test --workspace --all-targets

What we rejected. A single catch-all cargo test job. Splitting fmt / clippy / test / examples / msrv into separate jobs means each failure is reported independently and clearly, and the CI can tell you exactly which guarantee was broken.

4.5 A peek at the future: the sub-crate layout

The workspace’s members list is empty today, but the plan (PRD §3.1) is a folder tree like this:

flowchart TD
    accTitle: The planned sciencekit workspace layout
    accDescr: A sciencekit workspace root containing Cargo.toml, a crates folder with one sub-crate per algorithm area, and a docs folder.
    ROOT[sciencekit/]
    ROOT --> CARGO[Cargo.toml<br/>workspace root]
    ROOT --> CRATES[crates/]
    ROOT --> DOCS[docs/]
    CRATES --> UMBRELLA[sciencekit/<br/>umbrella crate re-exports]
    CRATES --> COMMON[sciencekit_common/<br/>base traits, types, errors]
    CRATES --> MATH[sciencekit_math/<br/>linear algebra, BLAS, SIMD, sparse]
    CRATES --> PREP[sciencekit_preprocessing/<br/>scalers, encoders]
    CRATES --> LINEAR[sciencekit_linear_model/<br/>regressions, logistic, SGD]
    CRATES --> MORE[... one crate per algorithm area]

    style ROOT fill:#e0f2fe,stroke:#0284c7,color:#0c4a6e
    style CRATES fill:#fef3c7,stroke:#d97706,color:#78350f
    style UMBRELLA fill:#f0fdf4,stroke:#16a34a,color:#14532d
    style COMMON fill:#f0fdf4,stroke:#16a34a,color:#14532d
    style MATH fill:#f0fdf4,stroke:#16a34a,color:#14532d

To see how a sub-crate actually registers into the workspace, here is the real sciencekit_math/Cargo.toml (from a later phase). Notice it is a member only once its own change adds it to the root’s members list, and it depends on sciencekit_common by path:

#![allow(unused)]
fn main() {
[package]
name = "sciencekit_math"
version = "0.1.0"
edition = "2024"
license = "Apache-2.0"
rust-version = "1.85"
description = "Computational kernels and the BLAS/LAPACK interface for sciencekit: pairwise distances, dense/sparse products, higher-order operations and memory-layout helpers."

[dependencies]
sciencekit_common = { path = "../sciencekit_common" }
faer = "0.24"
ndarray = { version = "0.17", features = ["rayon"] }
rayon = "1"
}

The path = "../sciencekit_common" is the workspace in action: an internal dependency resolved from a sibling folder, not from the internet.

5. A real walk-through, using the tests

Bootstrap (Phase 0.1) adds no algorithm code and therefore no *_tests.rs files — there is nothing to test yet. So the honest “executable example” here is the CI itself: what the pipeline actually runs on every change. Let us walk through the test job’s command as if we were the machine.

Step 1 — what cargo test --workspace --all-targets does

cargo test --workspace --all-targets

Run this in the terminal (it is what CI runs). Cargo:

1. Reads the root Cargo.toml and discovers the workspace and its members.
2. Compiles every crate in the workspace and their dependencies (the build).
3. Finds every #[test] function in every crate, builds them, and runs them.
4. Also builds and runs doctests — code examples written inside /// documentation comments. These are real, executable tests that prove the documented examples still work.
5. With --all-targets, it additionally compiles examples and benchmarks, catching compilation errors in those too, even though it only runs the tests.

At bootstrap the output is essentially “0 tests” — but the build succeeds, which proves the workspace itself is sound. That is the whole point of this phase: a green build with nothing inside is the green light to start adding crates.

Step 2 — how a gate blocks a change

Suppose a contributor (or an agent) pushes a change with a single formatting slip — say, an extra space. The fmt job runs:

cargo fmt --all -- --check

Because --check only verifies, it exits with a non-zero status code (a failure). GitHub Actions marks the job red. The pull request shows ❌ Formatting (cargo fmt). The change is blocked until the contributor runs cargo fmt --all (without --check) to actually fix the formatting, then pushes again.

This is the loop that keeps the codebase healthy forever: write → push → machine verifies every guarantee → fixes if any gate fails.

Step 3 — why the MSRV gate exists

Here is a scenario the msrv job is designed to catch. A developer’s local Rust is version 1.90. They use a brand-new API that only exists in 1.90. It compiles fine locally, and cargo test passes. But the project promises support for 1.85. If nothing checks, this promise silently breaks.

The msrv job installs exactly 1.85 and builds + tests again:

cargo build --workspace --all-targets
cargo test --workspace --all-targets

If that 1.90-only API was used, the build fails here — even though it passed on the developer’s machine. The promise is enforced, not assumed.

6. Look inside

The relationship between all the pieces we built:

flowchart TB
    accTitle: How the workspace, toolchain and CI gates fit together
    accDescr: The Cargo.toml workspace root ties together the toolchain file, the umbrella crate and the CI workflow. The crates subfolder grows one crate per change, and the CI pipeline runs fmt, clippy, tests, examples and an MSRV check as gates.
    Root["Cargo.toml<br/>(workspace root, members = [])"] --> Toolchain["rust-toolchain.toml<br/>channel = 1.85, edition 2024"]
    Root --> Umbrella["src/lib.rs<br/>#!forbid(unsafe_code)"]
    Root --> CI["ci.yml<br/>(GitHub Actions)"]

    subgraph Crates["crates/ (added one per change)"]
        Common["sciencekit_common"]
        Math["sciencekit_math"]
        More["... more later"]
    end

    Root -. grows into .-> Crates
    Common -. path dependency .-> Math

    subgraph Gates["CI gates"]
        Fmt["fmt"] --> Clippy["clippy -D warnings"]
        Clippy --> Test["test + doctests"]
        Test --> Ex["examples"]
        Ex --> Msrv["MSRV 1.85"]
    end
    CI --> Gates

    classDef root fill:#e0f2fe,stroke:#0284c7,color:#0c4a6e
    classDef crate fill:#fef3c7,stroke:#d97706,color:#78350f
    classDef gate fill:#f0fdf4,stroke:#16a34a,color:#14532d
    class Root,Toolchain,Umbrella,CI root
    class Common,Math,More crate
    class Fmt,Clippy,Test,Ex,Msrv gate

The workspace root points at the toolchain (what version to use), the umbrella lib.rs (the crate that will one day re-export everything), and the CI workflow. The CI runs the five gates in order. The crates/ folder starts empty and is filled one sub-crate per change, wired together by path dependencies.

7. Recap

• A Cargo workspace splits the library into small, isolated sciencekit_* sub-crates so each area compiles independently, keeps PRs small, and shares dependencies by path.
• Edition 2024 is a language compatibility mode; MSRV 1.85 is a compatibility promise. rust-toolchain.toml enforces 1.85 on our machines, while rust-version in Cargo.toml announces it to the world.
• CI is automated enforcement: the five gates (fmt, clippy -D warnings, test, examples, msrv) each catch a different failure class, and a red gate blocks the change until it is fixed.
• The naming convention (SK / sk_, no other abbreviations) is decided now and applied to every object in every later chapter.
• Bootstrap is the skeleton: it contains no algorithms, but every algorithm that follows stands on it.

Next chapter: with the workspace and CI in place, the first real content arrives — the shared contract vocabulary in sciencekit_common, where SKFloat, SKError, and the fundamental traits like SKEstimator begin.

Core contracts in sciencekit_common

Before any algorithm can crunch numbers, the whole library has to agree on a shared vocabulary: what a number is, what an error looks like, what a piece of data is, and what “fit” or “predict” even mean as a promise. That shared vocabulary lives in one crate, sciencekit_common. This chapter walks through the foundational traits and types every future algorithm will compile against — starting from a reader who has never seen Rust generics, traits, or even floating-point numbers.

1. The problem we were solving

Imagine a whole team of engineers each building a different algorithm: one builds a linear regression, another a decision tree, another a clustering method. If they each invent their own way to represent the data, their own error messages, and their own idea of what “fit” means, the library becomes a zoo — every algorithm speaks a different dialect, and nothing composes with anything else.

Rust makes this worse because of its philosophy of strict types. The compiler refuses to let two functions interoperate unless their types line up exactly. So the very first thing a library like this needs is not an algorithm at all — it is a contract: a set of shared building blocks that every future piece of code can agree on.

sciencekit_common is that contract. It is the crate that “freezes the vocabulary every future crate compiles against” (as its own lib.rs header says). It decides:

• what a valid floating-point number is (SKFloat),
• how failures are reported (SKError),
• how feature data and target data enter an algorithm (SKDataView, SKTargetView),
• what it means for an algorithm to fit, predict, or transform (SKSupervisedFit, SKUnsupervisedFit, SKRegressorPredictor, SKClassifierPredictor, SKFeatureTransformer).

Without this crate, sciencekit would be dozens of disconnected experiments. With it, every algorithm starts from a proven, shared foundation — like every room in a building sharing the same electrical standard and the same wall sockets.

Think of it as a plug-and-socket problem. Every appliance (algorithm) needs the same plug shape (the data-view contract), the same socket voltage (the scalar type), and the same way to tell you it broke (the error type). sciencekit_common manufactures those plugs, sockets, and fuses once, so every appliance can share a single wall.

2. The decisions — and the roads not taken

3. The concepts, taught from zero

This chapter leans on several Rust ideas and several numerical ideas. Let’s meet each before we look at real code.

Rust concept 1: what a trait is

A trait in Rust is like an interface in other languages: a promise that a type can do certain things. When we write

#![allow(unused)]
fn main() {
trait SKFloat: Float { /* ... */ }
}

we are saying: “anything that claims to be an SKFloat must also be a Float” — that is, it must support arithmetic like +, *, .abs(), and so on. Types then implement the trait:

#![allow(unused)]
fn main() {
impl SKFloat for f32 {}
impl SKFloat for f64 {}
}

This says “the types f32 and f64 promise to be SKFloat.”

Rust concept 2: generics — writing one function for many types

A generic function is written once but works on many types. In many languages you write template <typename T>; in Rust you write <T> or, in our case, <F: SKFloat>. The colon is a bound: “F can be any type, as long as it satisfies SKFloat.”

#![allow(unused)]
fn main() {
fn scale<F: SKFloat>(value: F, factor: F) -> F {
    value * factor
}
}

This function does not care whether F turns out to be f32 or f64 — it just knows that whatever F is, it supports * and can be returned. The caller chooses the type:

#![allow(unused)]
fn main() {
let x32 = scale(2.0_f32, 3.0_f32); // F is f32
let x64 = scale(2.0_f64, 3.0_f64); // F is f64
}

Both call the same function, but the compiler produces a specialized version for each concrete type. This is called monomorphization, and it is why generics cost nothing at runtime — there is no “boxing” or dynamic lookup; it is just ordinary statically typed code, duplicated per type.

Rust concept 3: sealed traits (why only f32/f64 can be SKFloat)

A sealed trait is a trait that outsiders cannot implement. We want SKFloat to be implemented only by f32 and f64 — not by some future third-party numeric type, and certainly not by integers. Rust has no built-in “seal” keyword, so we achieve sealing with a clever trick: make the trait require a private supertrait that only we can implement.

#![allow(unused)]
fn main() {
pub trait SKFloat: Float + Send + Sync + 'static + private::SKFloatSealed {}
}

The private::SKFloatSealed supertrait lives in a private module that is not exported. Because no outside crate can even name that module, no outside crate can write impl SKFloatSealed for MyType. And since implementing SKFloat requires implementing SKFloatSealed, only the two implementations we wrote (for f32 and f64) can ever exist. The sealing is enforced at compile time, by the module system itself.

This is subtle and powerful: the constraint is not checked at runtime — it cannot even be written by anyone outside the crate. That is exactly the safety we want.

Rust concept 4: enums — a value that is one of several shapes

An enum is a type whose value is exactly one of a fixed set of possibilities. In Rust each possibility can carry its own data. For example:

#![allow(unused)]
fn main() {
pub enum SKDataView<'a, F> {
    Dense(ArrayView2<'a, F>),   // a 2-D dense matrix view
    Sparse(CsMatView<'a, F>),   // a sparse matrix view
}
}

A SKDataView value is either a dense view or a sparse view — never both. When we later match on it, we write one branch per case and the compiler guarantees we handled every possibility. #[non_exhaustive] means “there may be more variants in the future, so don’t write a match without a wildcard _ arm.”

Rust concept 5: associated types (type Model)

Some traits promise not just methods but also a type that is chosen by the implementer. This is an associated type:

#![allow(unused)]
fn main() {
pub trait SKSupervisedFit<F: SKFloat> {
    type Model;   // the fitted model type, chosen by the implementer
    fn fit<'a, X, T>(&self, features: X, targets: T) -> Result<Self::Model, Self::Error>;
}
}

Whoever implements SKSupervisedFit must say what Model is. When a linear regression implements this trait, its Model is a fitted-linear-regression struct; when a decision tree does, its Model is a fitted-tree struct. The contract of fit is fixed, but the concrete fitted type is up to each algorithm. type Error is the same idea for errors: each algorithm may pick its own error type, as long as it can be built From the central SKError.

Rust concept 6: zero-copy views vs owned arrays

This is both a Rust and a performance concept. An owned array (Array2<f64>) owns its data in memory — it is the actual box of numbers. A view (ArrayView2<'a, f64>) is like a window that looks at somebody else’s numbers without owning them. The view is tiny (a few pointers and length fields), and it borrows the underlying buffer.

The 'a is a lifetime: it says “this view is only valid for as long as the data it borrows is alive.” Rust’s borrow checker guarantees you never hold a view to data that has been freed. The lib.rs crate even forbids unsafe_code at the crate level (#![forbid(unsafe_code)]), relying entirely on safe Rust and the borrow checker for memory safety.

Creating a view is free — no elements are copied. Copying a big matrix element-by-element is exactly the kind of waste this library was built to avoid.

Rust concept 7: the TryInto conversion seam

Rust’s standard library provides two conversion traits:

• Into<T> — an infallible conversion that cannot fail.
• TryInto<T> — a fallible conversion that returns a Result, so it can fail.

The key insight in sciencekit is that public APIs declare their inputs over TryInto, not over a fixed concrete type. So fit does not say “give me an Array2”; it says “give me anything that can TryInto a SKDataView”:

#![allow(unused)]
fn main() {
where
    X: TryInto<SKDataView<'a, F>, Error = SKError>,
}

Because f64: Into<T> is automatically also TryInto<T> (a blanket rule in the standard library), an infallible conversion is accepted for free. And a third-party library can integrate by writing its own impl TryFrom<TheirType> for SKDataView. This makes the library open to extension at its boundaries without ever changing the core.

Rust concept 8: error handling via Result + thiserror

Rust functions that can fail return a Result<T, E> — either Ok(value) on success or Err(error) on failure. The thiserror crate makes defining good error types easy: you write the enum, and #[derive(Error)] plus #[error("...")] attributes generate the Display and source implementations for you. That is why our error enum is so readable.

Numerical concept 1: what a floating-point scalar is

Numbers like 3.14 or -0.001 with fractional parts are floating-point numbers. A computer stores them in binary as sign + exponent + mantissa (scientific notation in base 2). They are approximate — not every decimal fraction can be represented exactly. The name “float” comes from the decimal point being able to “float” to any exponent.

Numerical concept 2: why f32 vs f64

• f32 is 32-bit single precision: smaller (4 bytes each), faster on some hardware, but stores about 7 significant decimal digits.
• f64 is 64-bit double precision: larger (8 bytes each), slower, but stores about 15–16 significant decimal digits — the standard for scientific computing.

Machine learning often defaults to f64 for correctness in training (where repeated arithmetic accumulates error) and can use f32 for memory/throughput when accuracy is good enough. Our sealed SKFloat allows both, so a single algorithm works with either — the caller picks, and the code is generic over both. Integers, by contrast, are not allowed here: continuous algorithms need fractional values and decimal precision, which integers fundamentally lack (spec scalar-typing).

Numerical concept 3: what a “view” over data means, and why zero-copy matters

A view over data is a lens, not a copy. If you have a 10-million-element matrix and you pass a view to a function, the function sees the same underlying numbers without the cost of duplicating 10 million elements. This is the heart of zero-copy: no bytes are moved; only a lightweight descriptor travels.

Why does zero-copy matter for performance? Copying memory is slow relative to arithmetic. In a pipeline of many stages (transform → scale → train), if every stage copied its input, a matrix would be duplicated over and over — each copy costing bandwidth and time. By borrowing views all the way through, the data stays in place and only descriptors move. On very large, possibly out-of-core datasets, avoiding copies is the difference between “fits in memory” and “swaps.”

Numerical concept 4: dense vs sparse (CSR) matrices

A dense matrix stores every cell, including all the zeros, in a regular grid. It is simple and fast when most cells are nonzero.

A sparse matrix stores only the nonzero entries plus enough bookkeeping to remember where they were. CSR stands for compressed sparse row: the matrix is stored row by row, keeping (a) the values, (b) which column each value lives in, and (c) pointers marking where each row begins. If 99% of a matrix is zeros (think user–item ratings, text, networks), storing everything densely wastes memory. Sparse storage packs the same data far more compactly.

A SKDataView can be either — an algorithm that works on both gets one uniform type to match on, instead of two code paths.

Numerical concept 5: what a “target” or “label” is in machine learning

In supervised learning you have features (the input, e.g. pixel values, measurements) and a target (what you want to predict — the answer). Three flavors matter here:

• Continuous: a real number you predict, e.g. house price → F, the model scalar (regression).
• Integer: a whole-number target, e.g. star ratings → i64 (elevated to F on demand).
• Nominal: a category, e.g. "cat", "dog", "bird" → strings (classification).

SKTargetView is the single type that holds any of these three. Unsupervised algorithms (e.g. clustering) have no targets at all — which is exactly why SKUnsupervisedFit::fit takes only features.

4. Each object, explained

SKFloat (the sealed floating-point bound)

What it is: The single trait that says “this type is a valid continuous number for sciencekit.” It is Float + Send + Sync + 'static + private::SKFloatSealed, implemented only for f32 and f64. Send + Sync + 'static mean the value can be moved across threads and lives forever (needed because algorithms may run in parallel).

Why it exists: Every generic algorithm needs a numeric bound, and having exactly one guarantees consistency (spec scalar-typing). Send + Sync + 'static ensure a value of type F is thread-transferable, which future parallel algorithms depend on.

What we rejected: A scatter of loose bounds on every function, and allowing third-party numeric types. Sealing locks the surface so only our two floats are valid.

private::SKFloatSealed (the private seal — non-public but essential)

What it is: A hidden supertrait, #[doc(hidden)], defined in a private module, with implementations only for f32 and f64.

Why it exists: It is the mechanism that seals SKFloat. Because the module is private, outsiders cannot name the trait, so they cannot implement it, and therefore cannot implement SKFloat. This is the compile-time guarantee behind the whole scalar contract.

What we rejected: Language-level sealed keywords don’t exist in Rust; this module trick is the idiomatic, documented way to achieve sealing.

SKError (the central error taxonomy)

What it is: A #[non_exhaustive] enum, derived with thiserror, with eight variants: ShapeMismatch, UnsupportedRepresentation, InvalidHyperparameter, ExecutionModeIncompatible, BatchBufferOversize, NotConverged, Io, and Conversion. Each carries precise data (e.g. expected/found shapes, iteration counts, byte counts) and a readable Display message.

Why it exists: One shared taxonomy means the same failure means the same error everywhere (spec error-model). The #[from] io::Error variant lets platform I/O errors convert in automatically; #[non_exhaustive] lets us grow the enum without breaking consumers.

What we rejected: Ad-hoc error types per algorithm and plain String errors — both lose precision or consistency. Convenience constructors like SKError::shape_mismatch_2d and SKError::not_converged exist to keep call sites tidy.

SKDataView (the zero-copy feature view)

What it is: A #[non_exhaustive] enum over the two feature representations: Dense(ArrayView2<'a, F>) and Sparse(CsMatView<'a, F>). Both are views — they borrow the caller’s data. It derives Debug, Clone, Copy.

Why it exists: One uniform type for “the input features,” whether dense or sparse, enables a single dispatch per operation (spec data-view-boundary). The .representation() method returns "dense"/"csr" for one-time dispatch and diagnostics; .as_dense() returns the dense view or rejects a sparse input before processing any element.

What we rejected: Passing owned arrays (forces copies) and separate functions per representation (duplication, no uniformity). #[non_exhaustive] keeps room for future variants.

SKTargetView (the zero-copy target view)

What it is: A #[non_exhaustive] enum over three target representations: Continuous(ArrayView1<'a, F>), Integer(ArrayView1<'a, i64>), and Nominal(&'a [&'a str]). Also derives Debug, Clone, Copy.

Why it exists: Supervised algorithms need a single type for “the answers” that covers regression (continuous), integer targets, and classification (nominal). .as_continuous() elevates continuous or integer targets to the model scalar F (exact while the value is representable in F) and rejects nominal with a clear suggestion.

What we rejected: Separate per-kind target arguments (forces the caller to know the kind up front) and forcing nominal text to be copied (we borrow the strings instead — spec data-view-boundary).

SKSupervisedFit (the supervised fit contract)

What it is: A trait with associated Model and Error types whose fit(&self, features, targets) takes both features and targets, each declared over the TryInto seam, and returns Result<Self::Model, Self::Error>.

Why it exists: It is the promise every supervised estimator keeps: given features and targets, produce a fitted model (spec estimator-contracts). Taking &self keeps the configured estimator unchanged and reusable.

What we rejected: Returning the learned state on the estimator itself (that would let you predict before fitting and prevent concurrent fits). The model is a separate type.

SKUnsupervisedFit (the unsupervised fit contract)

What it is: A trait with associated Model and Error types whose fit(&self, features) takes only features and returns Result<Self::Model, Self::Error>.

Why it exists: Clusterers and other unsupervised algorithms have no targets; making the signature lack targets makes it impossible to pass them (compile-time rejection). Taking &self keeps the estimator immutable and reusable, including concurrent fits.

What we rejected: Reusing one fit trait with an optional/ignored target parameter — that would let callers pass meaningless targets and lose the type-level guarantee.

SKRegressorPredictor (the regressor prediction contract)

What it is: A trait generic over F: SKFloat with an associated Error whose predict(&self, features) returns Result<ndarray::Array1<F>, Self::Error> over the TryInto seam.

Why it exists: Prediction lives on the fitted model, never on the configured estimator (spec estimator-contracts). That makes “predict before fit” a compile-time error. predict returns a dense Array1<F> — responses in the same scalar the model was fitted on, so an f32 pipeline never pays an f64 detour.

What we rejected: A single f64-only predict for every model. It forced an upcast on every f32 prediction and a cast back to chain the next stage — and it encoded class labels as floats (0.0/1.0), which is a lie of a type.

SKClassifierPredictor (the classifier prediction contract)

What it is: A trait generic over F: SKFloat with an associated Error whose predict_labels(&self, features) returns Result<ndarray::Array1<i64>, Self::Error> and whose predict_probabilities(&self, features) returns Result<ndarray::Array2<F>, Self::Error>, both over the TryInto seam.

Why it exists: Classifier outputs are labels, not floats. Labels are canonical i64 indices — independent of the feature scalar, so f32 and f64 runs agree exactly — while probabilities preserve the model scalar F.

What we rejected: Putting predict_labels on the estimator type. That would break the “cannot predict before fitting” guarantee and the shared model across threads.

SKFeatureTransformer (the feature-transformer contract)

What it is: A trait with associated Output and Error types whose transform(&self, features) returns Result<Self::Output, Self::Error> over the TryInto seam.

Why it exists: Transformers (e.g. scalers) need a typed Output so pipelines can check at compile time that one stage’s output feeds the next stage’s input (spec estimator-contracts). The associated type makes incompatible chaining fail early.

What we rejected: Returning an untyped or fixed output — that would push type compatibility checking to runtime.

Here is the whole family at a glance:

flowchart TD
    accTitle: The sciencekit_common vocabulary at a glance
    accDescr: Numbers, data views and errors feed the fit traits on the estimator and the prediction and transform traits on the fitted model.
    COMMON[sciencekit_common] --> NUM[numbers]
    COMMON --> DATA[data and errors]
    NUM --> FL[SKFloat sealed: f32, f64]
    FL --> SEAL[SKFloatSealed private seal]
    DATA --> ERR[SKError enum]
    DATA --> DV[SKDataView: Dense or Sparse]
    DATA --> TV[SKTargetView: Continuous, Integer or Nominal]
    TV --> SF[SKSupervisedFit: features plus targets]
    DV --> SF
    DV --> UF[SKUnsupervisedFit: features only]
    SF --> RP[SKRegressorPredictor on the model]
    SF --> CP[SKClassifierPredictor on the model]
    UF --> M[model type]
    DV --> RP
    DV --> CP
    RP --> FT[SKFeatureTransformer: transform]

    style COMMON fill:#e0f2fe,stroke:#0284c7,color:#0c4a6e
    style NUM fill:#ede9fe,stroke:#7c3aed,color:#4c1d95
    style DATA fill:#fef3c7,stroke:#d97706,color:#78350f
    style SF fill:#f0fdf4,stroke:#16a34a,color:#14532d
    style UF fill:#f0fdf4,stroke:#16a34a,color:#14532d
    style RP fill:#f0fdf4,stroke:#16a34a,color:#14532d
    style CP fill:#f0fdf4,stroke:#16a34a,color:#14532d
    style FT fill:#f0fdf4,stroke:#16a34a,color:#14532d

5. A real walk-through, using the tests

The most convincing way to learn is to read actual tests that passed. Here are three real tests from data_view_tests.rs and sk_float_tests.rs.

Test 1: borrowed_dense_enters_without_copying

This test proves the zero-copy promise: a dense matrix enters the library without copying.

#![allow(unused)]
fn main() {
#[test]
fn borrowed_dense_enters_without_copying() {
    let data = Array2::from_shape_vec((2, 2), vec![1.0_f64, 2.0, 3.0, 4.0]).unwrap();
    let view: SKDataView<'_, f64> = ArrayView2::from(&data).try_into().unwrap();
    match view {
        SKDataView::Dense(v) => {
            assert_eq!(v.as_slice_memory_order().unwrap(), &[1.0, 2.0, 3.0, 4.0]);
        }
        other => panic!("expected dense, got {other:?}"),
    }
}
}

Line by line:

1. let data = Array2::from_shape_vec((2, 2), vec![1.0_f64, 2.0, 3.0, 4.0]).unwrap(); — we build an owned 2×2 dense matrix (2 rows, 2 columns) holding the four values 1.0, 2.0, 3.0, 4.0. The _f64 suffix says the numbers are 64-bit floats. .unwrap() unwraps the Result; here it cannot fail.
2. let view: SKDataView<'_, f64> = ArrayView2::from(&data).try_into().unwrap(); — first we make a view of data (borrowing it, no copy), then we ask to convert it into a SKDataView. The '_ lifetime says “this view is valid for as long as data lives.” The conversion is the TryFrom<ArrayView2> impl we saw — it just wraps the view as SKDataView::Dense. No elements are copied.
3. match view { ... } — we check which variant we got. Since we handed it a dense view, we expect SKDataView::Dense(v).
4. assert_eq!(v.as_slice_memory_order().unwrap(), &[1.0, 2.0, 3.0, 4.0]); — we confirm the dense view exposes exactly the four original values, in the same order. The other => panic! arm catches the (impossible here) case where it came back sparse.

The lesson: the same buffer data is reused; the test’s very name asserts no copy was made.

Test 2: dense_only_consumer_rejects_sparse_precisely

This test proves that an algorithm supporting only dense data rejects a sparse input with a precise error, before touching any element.

#![allow(unused)]
fn main() {
#[test]
fn dense_only_consumer_rejects_sparse_precisely() {
    let csr = CsMat::new_csc((2, 2), vec![0, 0, 0], vec![], Vec::<f64>::new());
    let view: SKDataView<'_, f64> = csr.view().try_into().unwrap();
    let dense_result = view.as_dense();
    match dense_result {
        Err(SKError::UnsupportedRepresentation { representation, .. }) => {
            assert_eq!(representation, "csr");
        }
        other => panic!("expected unsupported-representation, got {other:?}"),
    }
}
}

Line by line:

1. CsMat::new_csc((2, 2), vec![0, 0, 0], vec![], Vec::<f64>::new()) — build a 2×2 sparse matrix with no stored values (an empty sparse matrix). Note the new_csc here really produces a sparse matrix in the sprs crate’s row-compressed form usable as CSR.
2. let view: SKDataView<'_, f64> = csr.view().try_into().unwrap(); — take a view of the sparse matrix and convert it through the same seam. This time it becomes SKDataView::Sparse(...).
3. let dense_result = view.as_dense(); — we simulate a dense-only algorithm calling the dispatch method. Since the view is sparse, as_dense() returns an Err.
4. The match checks that the error is exactly SKError::UnsupportedRepresentation with representation == "csr" — not a vague shape error. This matches the spec: “the operation fails before processing any element, with an error indicating the mismatch and the suggested conversion path.” The other => panic! arm catches any wrong error kind.

The lesson: representation errors are precise and distinct from shape errors, and they happen up front — never after a wasted pass over the data.

Test 3: native_floats_satisfy_the_contract

This test proves the sealed SKFloat works for both supported floats, through a generic function.

#![allow(unused)]
fn main() {
#[test]
fn native_floats_satisfy_the_contract() {
    let x32 = scale(2.0_f32, 3.0_f32);
    assert_eq!(x32, 6.0_f32);

    let x64 = scale(2.0_f64, 3.0_f64);
    assert_eq!(x64, 6.0_f64);
}
}

The helper scale is defined just above it:

#![allow(unused)]
fn main() {
fn scale<F: SKFloat>(value: F, factor: F) -> F {
    value * factor
}
}

Line by line:

1. fn scale<F: SKFloat>(value: F, factor: F) -> F — a generic function over F, bounded by SKFloat. It takes two values and returns one, all of type F.
2. value * factor — the body uses *. This compiles only because SKFloat: Float guarantees multiplication exists. The function never needs to know whether F is f32 or f64.
3. In the test, scale(2.0_f32, 3.0_f32) instantiates F = f32 and yields 6.0_f32; scale(2.0_f64, 3.0_f64) instantiates F = f64 and yields 6.0_f64. Both assert_eq! checks pass.

The companion test integers_do_not_satisfy_the_contract documents the other side: i32 does not implement SKFloat, so calling a continuous function with an integer is a compile-time error. And floats_are_thread_transferable asserts f32/f64 are Send + Sync + 'static, guaranteeing they can cross thread boundaries — the basis for future parallel algorithms.

6. Look inside

Here is how the fit/predict/transform contracts relate to the estimator, the model, and the data views:

flowchart TD
    accTitle: The fit, predict and transform contracts and their shared inputs
    accDescr: The configured estimator carries SKSupervisedFit and SKUnsupervisedFit, which take features and targets and produce a fitted model. The fitted model carries SKRegressorPredictor for continuous responses in the model scalar, SKClassifierPredictor for scalar-independent integer labels plus scalar-preserving probabilities, and SKFeatureTransformer. All of them consume SKDataView and SKTargetView, and all report errors through SKError.
    subgraph est[Configured estimator - immutable, reusable]
        SF[SKSupervisedFit] -->|fit features + targets| M1(Model type)
        UF[SKUnsupervisedFit] -->|fit features only| M2(Model type)
    end

    subgraph model[Fitted model - sole bearer of learned state, thread-shareable]
        RP[SKRegressorPredictor] -->|predict features| OUT[Array1 F]
        CP[SKClassifierPredictor] -->|predict labels| OUTL[Array1 i64]
        CP -->|predict probabilities| OUTP[Array2 F]
        FT[SKFeatureTransformer] -->|transform features| OUT2[Output type]
    end

    DV["SKDataView: Dense or Sparse"] --> SF
    DV --> UF
    DV --> RP
    DV --> CP
    DV --> FT
    TV["SKTargetView: Continuous, Integer or Nominal"] --> SF

    ER[SKError] --> SF
    ER --> UF
    ER --> RP
    ER --> CP
    ER --> FT

    style est fill:#e0f2fe,stroke:#0284c7,color:#0c4a6e
    style model fill:#f0fdf4,stroke:#16a34a,color:#14532d
    style DV fill:#fef3c7,stroke:#d97706,color:#78350f
    style TV fill:#fef3c7,stroke:#d97706,color:#78350f
    style ER fill:#ffe4e6,stroke:#e11d48,color:#7f1d1d

The diagram shows the two halves of the system: the estimator side (which fits) and the model side (which predicts/transforms). The fitted model is deliberately a separate type from the estimator, and both depend only on the shared data views and the central error.

7. Which float should I use?

Both f32 and f64 satisfy SKFloat, and a model fitted on one predicts in the same one — the scalar flows end to end with no implicit cast. The choice is a trade-off:

• f32 halves memory traffic and doubles throughput on SIMD and GPU kernels, which is what out-of-core and large-batch workloads want. Prefer it unless you need the extra precision.
• f64 keeps roughly 15–16 decimal digits of precision for ill-conditioned solves and long accumulations. Prefer it when numerical stability dominates.

A pipeline carries a single scalar: every stage agrees on F, so mixing an f32 stage into an f64 pipeline fails at compile time until an explicit cast stage bridges them.

Migration note

Before this change, predict returned Array1<f64> and continuous targets were ArrayView1<f64> regardless of the feature scalar. After it, regressor predictions and continuous targets use Array1<F>/ArrayView1<F>, and classifier labels are Array1<i64> instead of float-encoded indices. Call sites that assumed f64 update to the model scalar F; classifiers switch from comparing floats to comparing label indices.

8. Recap

• SKFloat is sealed — only f32 and f64 are valid continuous numbers, enforced at compile time by a private supertrait (private::SKFloatSealed).
• Predictions preserve the model scalar — regressors return Array1<F>, classifiers return Array1<i64> labels plus Array2<F> probabilities; continuous targets are SKTargetView::Continuous(ArrayView1<F>).
• SKError is one precise taxonomy — shape, unsupported representation, hyperparameter, execution mode, batch overflow, convergence, I/O, and conversion — built with thiserror, #[non_exhaustive], and automatic From conversion into per-algorithm errors.
• Data enters through zero-copy views — SKDataView (dense/sparse) and SKTargetView (continuous/integer/nominal) borrow the caller’s buffers rather than copying them, which is the heart of the performance promise.
• Public inputs are declared over TryInto — any type that can convert into a view is accepted, making the library open to third-party integration without breaking.
• Fit, predict, and transform are separate contracts — SKSupervisedFit and SKUnsupervisedFit split by supervision; SKRegressorPredictor, SKClassifierPredictor, and SKFeatureTransformer live on the model type, making “predict before fit” a compile-time error.

Next chapter: now that we have the shared vocabulary of numbers, errors, and data views, we can look at how an actual algorithm puts it to work — fitting a model and predicting on new data using these very contracts.

sciencekit_math: kernels, distances and the BLAS interface

Every machine-learning algorithm, from linear regression to clustering to a support-vector machine, is at its heart just arithmetic done very, very fast over big tables of numbers. Before sciencekit could train its first model, it needed a solid numeric engine — a crate that knows how to walk every element of a table, how to measure how far apart two rows are, how to multiply a sparse matrix by a dense one, and how to outsource the heaviest algebra to a tuned library. That engine is sciencekit_math. This chapter opens it up from the very beginning: first the Rust ideas it leans on, then the math, then the real objects and the real tests that prove they work.

1. The problem we were solving

Imagine you are teaching a machine to tell apart two types of flowers by measuring their petal lengths and widths. You store the measurements as a table: each row is one flower, each column is one measurement (a “feature”). To build most models you must do three kinds of work on such a table:

1. Transform every number. Normalize a column so every value sits between 0 and 1, or scale every value by some constant.
2. Summarize. Add up a column, or the whole table.
3. Measure similarity. “How far apart are these two flowers?” — the heart of clustering and nearest-neighbor methods.

That sounds simple, but doing it well is a performance minefield. A realistic dataset has millions of rows. If you loop over every element by hand, one by one, in the slowest way, even a tiny model takes ages. And the algorithms downstream don’t want to care how the arithmetic is done — they just want a number back.

The bigger problem is the heavy algebra. Models routinely need things like the inverse of a matrix, eigenvalues, or a singular-value decomposition — the “H-bombs” of linear algebra. Writing those from scratch is both slow and error-prone. Professional libraries have spent decades tuning them. But pulling in the traditional C/Fortran implementations (“system BLAS”) means the library will not even compile on a machine without a C toolchain — a heavy burden for a pure-Rust project.

So sciencekit_math had to solve two things at once:

• a fast, reusable set of small operations (the “kernels”) written the Rust way — with zero copying and parallelism — that downstream algorithms call like building blocks; and
• a single clean door (an interface) through which all the heavy algebra passes, so that today it uses a pure-Rust implementation and tomorrow, if a user enables it, a blazing system BLAS — without changing a single line of algorithm code.

The rest of this chapter is the story of how those two pieces were built and why they look the way they do.

2. The decisions — and the roads not taken

3. The concepts, taught from zero

This chapter needs a handful of Rust ideas and a handful of mathematical ones. Let’s meet them gently before touching any real code.

Rust concept 1: closures — functions you can hand around

In Rust you can write a small function “in place,” inline, right where you need it. Such an anonymous function is called a closure. You have already seen them in this book’s earlier chapters (|x| 2.0 * x + 1.0). The key superpower of a closure is that you can pass it as an argument to another function.

#![allow(unused)]
fn main() {
// `transform` is a parameter that is itself a function: it takes an F and returns an F.
pub fn sk_elementwise_transform<F: SKFloat>(
    input: &ArrayView<F, ndarray::Ix1>,
    transform: impl Fn(F) -> F + Sync + Send,
    parallelism: usize,
) -> Array<F, ndarray::Ix1>
}

Read the middle parameter type: impl Fn(F) -> F + Sync + Send. Fn is the Rust name for “this is a callable thing,” and (F) -> F says it takes one value of type F and returns one. So sk_elementwise_transform is a generic worker: it does not know what the operation is — the caller supplies the operation as a closure, and the worker applies it to every element. This is the essence of a higher-order function: a function that takes or returns other functions. It separates “what computation” (the closure, chosen by the caller) from “how to sweep the array” (the kernel, written once). The + Sync + Send bits are promises that the closure can be shared across threads — which we need when we parallelize.

Rust concept 2: azip! and par_azip! — elementwise iteration done right

The whole library refuses to write for i in 0..n manual index loops for elementwise work. Instead it uses the ndarray macros azip! and par_azip!, which iterate element-by-element across one or more arrays in lockstep, element 0 with element 0, element 1 with element 1, and so on. Look at the body of the transform:

#![allow(unused)]
fn main() {
azip!((out in &mut output, a in input) { *out = transform(*a); });
}

This says: “walk input and output together; for each pair, put transform(*a) into *out.” The macro handles the striding and indexing for you. par_azip! does the same thing in parallel — it splits the work across a thread pool. Because both macros share the exact same body syntax, you get parallelism almost for free: the only change between the sequential and parallel versions is the macro name.

Rust concept 3: owned arrays vs zero-copy views

An owned array, Array<F, Ix2>, owns its numbers in memory — it is the sole holder and can resize, mutate, or be dropped. A view, ArrayView2<F>, does not own data; it is a “window” that borrows someone else’s owned array, remembering the shape and how to reach each element. Views are the library’s secret to speed: a function can take a view of your existing buffer, do work, and hand you an answer — without ever copying your data. Copying a matrix of a million rows is expensive; borrowing a view of it is essentially free. That is why every public kernel here takes &ArrayView2<F> (a borrow of a view) and only returns a fresh owned Array2<F> when it genuinely needs a new result.

Rust concept 4: traits and a trait-object-free generic backend

A trait is a promise that a type can do certain things (an interface). Here the central trait is SKMathBackend:

#![allow(unused)]
fn main() {
pub trait SKMathBackend<F: SKFloat>: Send + Sync {
    fn gemm(&self, a: ArrayView2<F>, b: ArrayView2<F>, alpha: F, parallelism: usize) -> Array2<F>;
    fn svd(&self, a: ArrayView2<F>) -> Result<SKSingularValueDecomposition<F>, SKError>;
    // ... more: qr, cholesky, solve, eigh, lu, ...
}
}

Any type that impls this trait is a “math backend”: it promises to perform GEMM (matrix multiplication), SVD, QR, and so on. The trick here is that the trait is generic over F and the concrete backend type is known at compile time — the caller writes let backend = SKFaerBackend::new(); and then calls backend.gemm(...). Because the concrete type is fixed, the compiler generates a direct, specialized call. This is the opposite of a trait object (a Box<dyn SKMathBackend>), which would defer the choice to runtime and pay for a dynamic dispatch on every call. For a performance crate, avoiding that indirection — “trait-object-free” — is deliberate: heavy algebra runs on hot paths where an extra virtual call per operation would add up.

Numerical concept 1: what a “kernel” is

In numerical libraries a kernel is a small, reusable operation applied over an array or matrix. We meet three families:

• Elementwise transform — a function applied to every element independently: x → 2x + 1. Input and output have the same shape.
• Reduction — collapse many values into fewer (often one): a sum turns a column into a single number.
• In-place transform — mutate the existing buffer rather than allocate a new one, saving memory.

The sciencekit_math kernels are the “LEGO bricks” that downstream algorithms snap together.

Numerical concept 2: what a distance matrix is, and squared-Euclidean math

A distance matrix is a table where entry [i, j] is “how far apart are row i of the first set and row j of the second set.” Each row of a matrix is a point in multi-dimensional space (here, one point per flower).

The Euclidean distance between two points $a$ and $b$ is the length of the straight line between them: $\sqrt{\sum (a_k - b_k)^2}$. Its square is $\sum (a_k - b_k)^2$. Computing that literally means subtracting each coordinate pair. But algebra gives us a shortcut by expanding the square:

$$ |a - b|^2 = \sum (a_k - b_k)^2 = \sum (a_k^2 - 2 a_k b_k + b_k^2) = |a|^2 - 2(a \cdot b) + |b|^2 $$

where $|a|^2 = \sum a_k^2$ (the squared “length” of $a$) and $a \cdot b = \sum a_k b_k$ is the dot product of $a$ and $b$. Why bother? Because we can precompute $|a|^2$ for every row once, and $|b|^2$ for every row of the other set once, and then each cell only needs one dot product — a tight, SIMD-friendly loop over contiguous memory. This is the norm-expansion formulation the spec mandates, and it is the whole reason sciencekit_math computes distances this way.

Numerical concept 3: cosine and Manhattan distances

Cosine distance measures the angle between two points rather than their distance. The cosine of the angle between $a$ and $b$ is $\frac{a \cdot b}{|a| , |b|}$. If the vectors point the same way the cosine is $1$; if they are perpendicular (orthogonal) it is $0$. Distance is then defined as $1 - \text{cosine}$, so identical directions score $0$ and orthogonal directions score $1$. Nice property: cosine distance is scale-invariant — $[1,2]$ and $[10,20]$ point the same way, so their cosine distance is $0$.

Manhattan distance (also called L1) is the distance you would walk along a city grid — only left/right and up/down, no diagonals: $\sum |a_k - b_k|$. It is the sum of the absolute coordinate differences.

Numerical concept 4: sparse data and the CSR format

A matrix is sparse when most of its entries are zero — think of text classification, where each row is a document and each column says “how many times word X appeared.” Storing all those zeros as real numbers wastes memory and time. The CSR (Compressed Sparse Row) format stores only the non-zero values, in three parallel lists:

• the non-zero values themselves,
• the column each value belongs to,
• a pointer array that says where each row begins.

To multiply a CSR matrix by a dense matrix, you only touch the non-zero entries — no need to ever build a big dense matrix full of zeros. sciencekit_math consumes CSR through a view (CsMatView) so even the sparse input is never copied.

Numerical concept 5: C vs F memory contiguity

A matrix is stored in memory as a long list of numbers. Row-major (C) order lays out each row back-to-back; column-major (F) order lays out each column back-to-back. Most languages default to row-major (the “C” convention); Fortran-derived libraries default to column-major (“F”). For a hot loop it matters a lot which layout a buffer has, because CPUs love reading memory that is sequential. A contiguous buffer (all elements one after another) is much faster to sweep with SIMD than a strided one. So the layout helpers detect which case we are in and, when a kernel needs contiguous memory, produce a contiguous copy only when necessary — borrowing when already contiguous (zero copy) and owning a fresh copy only when strided.

Numerical concept 6: what BLAS and LAPACK are

BLAS (Basic Linear Algebra Subprograms) is a de-facto standard set of routines for vector and matrix arithmetic — the crown jewel being GEMM, the general matrix multiplication $C = \alpha \cdot A \cdot B$. LAPACK builds on BLAS to provide the higher-level decompositions: SVD, QR, Cholesky, LU, eigenvalues, least squares. For decades these have been hand-tuned in C/Fortran (OpenBLAS, Intel MKL, and so on) to be staggeringly fast. The catch: they need a native C toolchain to compile. sciencekit_math solves this by defining its own SKMathBackend interface and letting the default implementation be faer, a modern pure-Rust BLAS/LAPACK — fast, and compiles anywhere Rust does.

4. Each object, explained

SKMathBackend (the dense algebra trait)

What it is: The single host-centric interface for all heavy dense algebra. It is a trait generic over SKFloat, whose methods take zero-copy ArrayView2 inputs and return owned Array2 results. The surface spans gemm, svd, qr, cholesky, solve_triangular, solve, eigh, lu, slogdet, pinv, inv, lstsq, and norm/vector_norm. It also has default implementations: for example slogdet is derived from LU, pinv and lstsq are composed from an SVD, and norm reduces host-side — so a backend can implement only the primitives and still get everything else for free.

Why it exists: Downstream algorithm crates must not care which library does the heavy algebra. They call methods on the trait; swapping faer for ndarray-linalg (under the blas-backend feature) is a config change, not a rewrite (spec blas-interface). The zero-copy views on input keep the promise of not copying the caller’s data.

What we rejected: Letting a concrete backend type leak into the algorithm surface, and using a trait object (Box<dyn ...>) — the latter would add runtime dispatch on hot numerical paths, which the performance mandate rules out.

SKFaerBackend (the pure-Rust default backend)

What it is: A tiny, Debug + Clone + Copy + Default unit struct — it holds no data — created with SKFaerBackend::new(). It implements SKMathBackend<f64> and SKMathBackend<f32> by delegating GEMM to faer’s blocked matmul (honouring the plan’s parallelism) and the decompositions to faer’s high-level SVD/QR/Cholesky/LU/eigen solvers. A sibling SKMatrixMultiplyBackend routes GEMM to matrixmultiply and the decompositions to faer.

Why it exists: The default must work with no C toolchain. faer compiles on Rust 1.85 and is competitive, so the default build links no system BLAS at all.

What we rejected: Making a system BLAS the default, and oxiblas (fails to compile on the project’s Rust version). The unsafe here is contained in the SIMD/GEMM glue only.

sk_elementwise_transform (the elementwise kernel)

What it is: A free function that applies a closure x → transform(x) to every element of a 1-D view and returns a new owned array. It reads the resolved parallelism: 1 runs the sequential azip! form (zero dispatch overhead); > 1 runs par_azip! across the compute pool.

Why it exists: It is the most basic LEGO brick — every scaling, normalization or feature transform reuses it (spec higher-order-kernels). Iteration is delegated to azip!, never a manual index loop.

What we rejected: Hand-written index loops, and always-owning a fresh buffer when the caller might prefer an in-place variant (that is sk_scale_in_place’s job instead).

sk_binary_combine (the pairwise elementwise kernel)

What it is: A free function that pairs two equal-length views element by element, combine(left_k, right_k), and returns a new owned array. It asserts the lengths match and, like the transform, honours parallelism (1 → zip_mut_with, > 1 → par_azip!).

Why it exists: Many operations need two arrays at once — adding a bias term, multiplying two feature sets, comparing two signals. It is the binary sibling of the transform.

What we rejected: Requiring a pre-allocated output buffer from the caller; the function owns its result, keeping call sites simple.

sk_axis_sum (the reduction kernel)

What it is: A free function that sums a 2-D array along one axis into a 1-D result. axis = 0 collapses rows → per-column sums (length = number of columns); axis = 1 collapses columns → per-row sums (length = rows). It honours parallelism. The parallel axis-0 path splits rows into chunks, reduces each chunk into its own accumulator, then combines the partials — deliberately avoiding any shared mutable accumulator, so there is no data race.

Why it exists: Summing along an axis is everywhere (totals, statistics, normalization denominators). The race-free parallel design is the careful part (spec higher-order-kernels).

What we rejected: Writing directly into a single shared accumulator from multiple threads (that would race), and manual index loops for the reductions.

sk_scale_in_place (the in-place kernel)

What it is: A free function that multiplies every element of a 2-D array by a scalar factor in place, returning nothing. It mutates the caller’s buffer, allocating no new array (spec higher-order-kernels).

Why it exists: When you only need to scale existing data, copying the whole array is wasteful; in-place is the zero-allocation path.

What we rejected: Returning a new array for scaling (that is sk_elementwise_transform’s mode); this function’s whole point is to avoid the copy.

SKMemoryLayout, sk_memory_layout, sk_force_contiguous (the layout helpers)

What they are: SKMemoryLayout is a three-variant enum — CContiguous (row-major), FContiguous (column-major), Strided (neither). sk_memory_layout detects which case a view is in; sk_is_c_contiguous / sk_is_f_contiguous ask one yes/no question; and sk_force_contiguous returns a CowArray: a zero-copy borrow when the input is already contiguous, or a fresh owned contiguous copy when it is strided.

Why they exist: Hot SIMD kernels want contiguous memory. Detecting layout lets code feed fast paths with contiguous buffers and only pay for a copy when it truly has to.

What we rejected: Always copying (defeats zero-copy), and assuming a layout that the data does not have (would silently run slow or wrong).

The distance functions: sk_squared_euclidean_distance_matrix, sk_euclidean_distance_matrix, sk_manhattan_distance_matrix, sk_cosine_distance_matrix

What they are: Four free functions that take two &ArrayView2<F> feature matrices and return an owned Array2<F> distance matrix. Squared-Euclidean uses the norm-expansion $|a|^2 - 2a \cdot b + |b|^2$, with row norms precomputed and each cell needing one SIMD dot product; Euclidean is its element-wise square root; Manhattan is $\sum |a_k - b_k|$; cosine is $1 - \frac{a \cdot b}{|a| , |b|}$, with a guard that maps zero-norm rows to distance $1$ (matching scikit-learn).

Why they exist: Nearest-neighbor search, clustering, SVM kernels and imputation all need distance matrices. The norm-expansion formulation is the performance heart (spec pairwise-distances).

What we rejected: The naive per-pair subtract loop $\sum (a_k - b_k)^2$ for squared Euclidean (slower, not SIMD-friendly), and dividing by zero in cosine without a guard.

sk_csr_dense_product and sk_sparse_product (the sparse products)

What they are: sk_csr_dense_product multiplies a CsMatView by a dense ArrayView2, walking each non-zero entry of the CSR operand once and accumulating value * dense_row into the output row via azip! — the sparse operand is never densified. sk_sparse_product multiplies two CsMatView<f64> and delegates to sprs’ sparse matrix-matrix product (SMMP), returning a sparse CsMat<f64>.

Why they exist: Sparse workloads (Lasso, linear SVMs, text classification) would blow up memory if densified (spec dense-sparse-products). Both take zero-copy views.

What we rejected: Densifying the sparse operand before multiplying, and copying the sparse input.

5. A real walk-through, using the tests

The most convincing way to learn is to read actual tests that passed. Here are three real tests: one from kernels_tests.rs, one from pairwise_tests.rs, and one that proves the parallel and sequential paths agree.

Test 1: axis_one_sum_matches_row_totals (from kernels_tests.rs)

This test checks that sk_axis_sum with axis = 1 sums each row correctly.

#![allow(unused)]
fn main() {
#[test]
fn axis_one_sum_matches_row_totals() {
    let input: Array2<f64> = array![[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]];
    let sums = sk_axis_sum(&input.view(), 1, 1);
    assert_eq!(sums, array![6.0, 15.0, 24.0]);
}
}

Line by line:

1. let input: Array2<f64> = array![[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]; — we build an owned 3×3 matrix of f64. The array! macro lays out three rows, each with three numbers, so input has 3 rows and 3 columns.
2. let sums = sk_axis_sum(&input.view(), 1, 1); — we call the kernel. The &input.view() hands the function a view (a borrow, no copy). The first 1 is the axis: we are summing along axis 1, which collapses columns and keeps rows. The second 1 is the parallelism: run sequentially.
3. assert_eq!(sums, array![6.0, 15.0, 24.0]); — we compare to the hand-computed row totals. Row [1,2,3] sums to 6, row [4,5,6] to 15, row [7,8,9] to 24. The kernel’s sequential path iterates rows and folds each one, so the output is exactly this.

Test 2: squared_euclidean_matches_direct_computation (from pairwise_tests.rs)

This test proves the norm-expansion squared-Euclidean kernel against hand-computed values.

#![allow(unused)]
fn main() {
#[test]
fn squared_euclidean_matches_direct_computation() {
    let left: Array2<f64> = array![[0.0, 0.0], [3.0, 4.0]];
    let right: Array2<f64> = array![[1.0, 1.0]];
    let distance = sk_squared_euclidean_distance_matrix(&left.view(), &right.view());
    // dist([0,0],[1,1]) = 2 ; dist([3,4],[1,1]) = 4 + 9 = 13.
    assert!((distance[[0, 0]] - 2.0).abs() < 1e-12);
    assert!((distance[[1, 0]] - 13.0).abs() < 1e-12);
}
}

Line by line:

1. let left: Array2<f64> = array![[0.0, 0.0], [3.0, 4.0]]; — the first feature set has two points: (0,0) and (3,4). let right: Array2<f64> = array![[1.0, 1.0]]; — the second set has one point (1,1). So the distance matrix will be 2×1.
2. let distance = sk_squared_euclidean_distance_matrix(&left.view(), &right.view()); — we hand both sets in as views. Inside, the kernel computes $|a_{\text{row}}|^2$ for both rows ($0$ and $25$), $|b_{\text{row}}|^2$ for the single right row ($2$), then fills the $2 \times 1$ matrix with $|a|^2 + |b|^2 - 2(a \cdot b)$ per cell.
3. assert!((distance[[0, 0]] - 2.0).abs() < 1e-12); — the distance from (0,0) to (1,1) is $0 + 2 - 2(0 \cdot 1 + 0 \cdot 1) = 2$, matching the straight-line square $(1-0)^2 + (1-0)^2 = 2$. The .abs() < 1e-12 is a floating-point tolerance: we expect exactly 2, allowing tiny rounding.
4. assert!((distance[[1, 0]] - 13.0).abs() < 1e-12); — from (3,4) to (1,1) the squared distance is (1−3)² + (1−4)² = 4 + 9 = 13, which the norm-expansion also produces.

Test 3: elementwise_parallel_matches_sequential_on_large_array (from kernels_tests.rs)

This test checks that the parallel kernel agrees with the sequential one, which is the whole point of the parallelism switch.

#![allow(unused)]
fn main() {
#[test]
fn elementwise_parallel_matches_sequential_on_large_array() {
    let n = 1 << 16;
    let input: Array1<f64> = Array1::from_shape_fn(n, |i| (i as f64) / 7.0);
    let f = |x: f64| (x * 3.0 + 2.0).sin();
    let sequential = sk_elementwise_transform(&input.view(), f, 1);
    let parallel = sk_elementwise_transform(&input.view(), f, 8);
    assert_close(&sequential, &parallel);
}
}

Line by line:

1. let n = 1 << 16; — n is 2¹⁶ = 65536, a large array so that parallelism actually matters.
2. let input: Array1<f64> = Array1::from_shape_fn(n, |i| (i as f64) / 7.0); — we build a 1-D array of 65536 elements; element i is i / 7.0. from_shape_fn fills the shape by calling the closure for each index.
3. let f = |x: f64| (x * 3.0 + 2.0).sin(); — the transform: sin(3x + 2). Note the same closure is passed to both calls; this is the higher-order design in action.
4. let sequential = sk_elementwise_transform(&input.view(), f, 1); — run it sequentially (parallelism = 1).
5. let parallel = sk_elementwise_transform(&input.view(), f, 8); — run the same transform with parallelism = 8, which dispatches par_azip! across the pool.
6. assert_close(&sequential, &parallel); — the test’s own helper checks the two agree within a relative tolerance. Why tolerance instead of exact equality? Parallel reordering can shift the last bits of floating-point results, so we allow tiny differences. This is the contract that parallel kernels stay correct.

6. Look inside

Here is how the crate’s modules hang together:

flowchart LR
    accTitle: Modules of the sciencekit_math crate and their key functions
    accDescr: sciencekit_math groups kernels, layout, pairwise distances, sparse operations and a backend. Each module exposes free functions, and the backend module exposes the SKMathBackend trait implemented by faer, matrix multiply and ndarray-linalg.
    MATH[sciencekit_math] --> K[kernels]
    MATH --> L[layout]
    MATH --> PW[pairwise]
    MATH --> SO[sparse_ops]
    MATH --> BE[backend]

    K --> K1[sk_elementwise_transform]
    K --> K2[sk_binary_combine]
    K --> K3[sk_axis_sum]
    K --> K4[sk_scale_in_place]
    K --> AZ["azip! / par_azip! (higher-order iteration)"]

    L --> L1["SKMemoryLayout: C, F, or Strided"]
    L --> L2["sk_force_contiguous returns a CowArray"]

    PW --> P1[sk_squared_euclidean_distance_matrix]
    PW --> P2[sk_euclidean_distance_matrix]
    PW --> P3[sk_manhattan_distance_matrix]
    PW --> P4[sk_cosine_distance_matrix]
    P1 --> SIMD["wide f64x4 SIMD dot (simd_dot.rs)"]

    SO --> S1["sk_csr_dense_product (never densifies)"]
    SO --> S2["sk_sparse_product uses sprs SMMP"]

    BE --> TRAIT[SKMathBackend trait]
    TRAIT --> F[SKFaerBackend - pure Rust, default]
    TRAIT --> G[SKMatrixMultiplyBackend]
    TRAIT --> H["ndarray-linalg (opt-in blas-backend)"]

    style MATH fill:#e0f2fe,stroke:#0284c7,color:#0c4a6e
    style K fill:#fef3c7,stroke:#d97706,color:#78350f
    style L fill:#fef3c7,stroke:#d97706,color:#78350f
    style PW fill:#f0fdf4,stroke:#16a34a,color:#14532d
    style SO fill:#f0fdf4,stroke:#16a34a,color:#14532d
    style BE fill:#ffe4e6,stroke:#e11d48,color:#7f1d1d
    style TRAIT fill:#ffe4e6,stroke:#e11d48,color:#7f1d1d

The diagram shows sciencekit_math as the numeric substrate: the four module families (kernels, layout, pairwise, sparse_ops) supply the reusable operations, all built on azip!/par_azip! and zero-copy views, while backend supplies the single heavy-algebra door that defaults to pure-Rust faer and can be swapped for a system BLAS via a feature flag.

7. Recap

• Kernels separate what from how. sk_elementwise_transform and sk_binary_combine take a closure (the operation) and sweep memory with azip!/par_azip! — so switching to parallel is a one-word change, never a rewrite.
• Reductions are race-free by design. sk_axis_sum accumulates per-chunk partials in parallel rather than writing into one shared accumulator, so parallel results match the sequential ones within a tiny floating-point tolerance.
• Distances use math to win. Squared-Euclidean is computed as $|a|^2 - 2a \cdot b + |b|^2$, reusing a SIMD dot product per cell; cosine and Manhattan have their own clean formulas, all over zero-copy row views.
• Sparse stays sparse. sk_csr_dense_product walks only the non-zero entries and never densifies, while sk_sparse_product delegates to sprs’ SMMP.
• Heavy algebra hides behind one trait. SKMathBackend with the pure-Rust SKFaerBackend default, plus an opt-in blas-backend feature for a system BLAS — the default needs no C toolchain, and algorithms never change when the backend does.

Next chapter: armed with kernels, distances, sparse products and a linear-algebra door, we are ready to watch the first real estimator put them to work — fitting a model and making predictions on new data.

Execution modes, resolution and observability

A machine-learning library cannot assume one-size-fits-all execution: the same algorithm must run instantly on a small in-memory dataset, stream through one too large to fit in RAM, or shuffle bytes off a memory-mapped file. And when it does run, an engineer needs to see what happened — which mode was chosen, how long it took, whether it failed. This chapter explains how sciencekit turns a loose intent (“just run it the best way”) into a concrete plan (a specific mode, thread count, and batch size), how that plan drives a double-buffered streaming pipeline, and how every operation emits a trace span you can export to OpenTelemetry.

1. The problem we were solving

Imagine you build a library that does math on datasets. Some datasets are tiny — a few hundred rows, the kind you could hold in a spreadsheet. Others are enormous — a terabyte of sensor readings that will never, ever fit in your computer’s memory at once.

If you write your algorithm assuming one of these, you break the other. Assume “everything fits in memory” and the big dataset crashes or swaps painfully. Assume “always stream” and the tiny dataset pays an absurd overhead of threading and batching for no reason.

The naive solution is to make the caller decide: pass a flag that says “in-memory” or “streaming”. That works, but it is a burden on the user — and, worse, a source of bugs. The user might pick “in-memory” for data that is actually too large, or pick “streaming” for an algorithm that fundamentally needs random access (jumping to arbitrary positions), which streaming cannot provide.

The library should be smarter than that. It should let the user say “just pick the right one” (an automatic intent), and then decide for itself, based on how much memory is available, how big the dataset is, and what kind of access the algorithm needs. That is the first half of this chapter: execution mode resolution.

The second half is observability. When something runs — especially something slow or concurrent — you want to know what happened. Did fit choose in-memory or streaming? How many milliseconds did it take? Did it fail? Answering these questions requires wrapping every operation in a “span” — a named, timed box of structured metadata — that a tracing system can collect and export. This must be cheap by default and opt-in for export, so the default library never sends telemetry traffic over the network.

Both halves share one philosophy, which the previous chapter called the library’s core strength: push mistakes out of the running program and into the compiler and the types. We will see that philosophy again and again — in enums that make invalid states unrepresentable, in a deterministic resolver that cannot depend on hidden global state, and in a builder pattern that makes an unconfigured model impossible to construct.

2. The decisions — and the roads not taken

3. The concepts, taught from zero

Rust concept 1: enums with data-carrying variants

An enum is a type whose value is exactly one of a fixed set of possibilities. In Rust each possibility (a variant) can carry its own data. We met this in chapter 2 with SKDataView. Here it powers the execution intent:

#![allow(unused)]
fn main() {
pub enum SKExecutionMode {
    Automatic,
    InProcessSynchronous,
    InProcessAsynchronous,
    OutOfCoreStreaming,
    OutOfCoreMemoryMapped,
}
}

A value of this type is one of those five — never two, never none. The compiler knows the full list, so when we match on it we must handle every variant, and the compiler refuses to let us forget one. #[non_exhaustive] above the enum tells the compiler: “outsiders should still add a wildcard arm, because we may add variants in the future.” Inside our own crate we are exhaustive; outside, consumers write a _ => catch-all to stay future-proof.

Rust concept 2: the builder pattern — mandatory here

Rust structs often hold many optional settings. The builder pattern solves “too many optional parameters” by giving you a chainable object that you configure step by step, then finally ask to produce the real value.

In sciencekit this is mandatory: you cannot call a constructor directly. Every public estimator exposes a builder, and construction goes through build() (spec base-builders). Two shared pieces make this uniform across every algorithm:

• SKBuilderState — the reusable “storage drawer” that remembers the execution intent and knows how to resolve it into a plan.
• SKBuilder — the trait (an interface) every builder implements, guaranteeing the same two methods: execution_mode(...) and build().

The reason this is mandatory is safety-by-construction: a model only ever exists after build() has validated its configuration and returned Result. You cannot accidentally use an unbuilt, half-configured model.

Rust concept 3: fluent setters and the Automatic default

A builder method that returns &mut Self lets you chain calls. The key default is that the execution intent starts as Automatic:

#![allow(unused)]
fn main() {
pub fn execution_mode<E>(&mut self, mode: E) -> &mut Self
where
    E: TryInto<SKExecutionMode>,
}

Notice the argument is generic over E: TryInto<SKExecutionMode>. We learned TryInto in chapter 2. The point here: you can pass an SKExecutionMode directly (it converts infallibly), and a mis-typed string cannot compile, because a string cannot TryInto an SKExecutionMode. The &mut self return means you can chain: .execution_mode(...) then .build().

Because SKBuilderState::new() starts at Automatic, a consumer who never calls execution_mode still gets automatic behavior — the requirement “a consumer who does not set it still gets automatic execution.”

Rust concept 4: Result and error propagation

Functions that can fail return Result<T, E> — either Ok(value) or Err(error). The resolver returns Result<SKExecutionPlan, SKError>, and the builder’s build() returns Result<Model, SKError>. SKError is the central taxonomy we met in chapter 2; two variants matter here:

#![allow(unused)]
fn main() {
#[error("execution mode `{mode}` is incompatible with the declared access pattern `{pattern}`")]
ExecutionModeIncompatible { mode: &'static str, pattern: &'static str },

#[error("streaming batch buffer overflow: {required_bytes} bytes needed (2 × {batch_bytes} per batch) exceeds {available_bytes} bytes available")]
BatchBufferOversize { batch_bytes: u64, required_bytes: u64, available_bytes: u64 },
}

The #[error("...")] comes from thiserror, which generates the human-readable message. The fields give the caller structured information — not just “it failed”, but exactly which mode, which pattern, how many bytes were needed. The ? operator propagates a Result up the call stack automatically, but here the resolver builds its errors by hand with Err(...) when it detects an incompatibility or an oversized buffer.

Rust concept 5: threads, Arc, Mutex, and Condvar (the streaming driver)

The streaming driver runs two threads that share a buffer. Rust’s answer to shared mutable state is precise:

• Arc<T> (atomic reference count) is a smart pointer that many threads can hold a clone of; when the last clone drops, the value is freed. It lets the producer thread and the consumer thread both “hold onto” the same shared state.
• Mutex<T> wraps a value so only one thread can read/write it at a time (mutual exclusion). You lock it to get a guard, work, then drop the guard to release it.
• Condvar (condition variable) is a way for a thread to wait until a condition becomes true and be notified when it changes. The producer waits when the buffer is full; the consumer waits when the buffer is empty.

You will see these three working together in the walk-through: the buffer is Mutex-protected state plus a Condvar for signaling, all shared through an Arc.

Numerical concept 1: what “in-memory vs out-of-core” means

• In-memory execution: the whole dataset sits in RAM, so the algorithm can touch any element at any time (random access) and pass over it repeatedly. Fast, but bounded by memory.
• Out-of-core execution: the dataset is bigger than memory, so it lives on disk and is processed in pieces. Two flavors: streaming (read a sequential batch, compute, move on) and memory-mapped (the OS lazily pages parts of the file into memory as you access them, giving random access to a file larger than RAM).

Numerical concept 2: what an access pattern is

An algorithm declares how it wants to touch data:

• Sequential: read a batch, move to the next, never look back — perfect for streaming.
• RandomAccess: jump to element i in O(1) — only possible in-memory or memory-mapped.
• Iterative: pass over the whole dataset repeatedly (e.g. gradient descent epochs) — fits in-memory or streaming, but not a single forward scan.

This is the crux of the incompatibility rule: you cannot stream an algorithm that needs random access, because a streaming pipeline only hands you batches in order.

Numerical concept 3: double buffering

In streaming, the reader (I/O) and the compute (math) are two different speeds. If the reader waits for compute to finish before reading the next batch, the I/O sits idle. If compute waits for the reader, the CPU sits idle.

Double buffering keeps both busy: while compute works on batch k, the reader fetches batch k+1 into a second buffer. When compute finishes batch k, the already-fetched k+1 is ready — no waiting. The buffer_depth in the plan (double buffering = 1) is the number of prefetched batches ahead. The guard in the resolver ensures two resident batches fit in memory: if one batch is $B$ bytes, double buffering needs $2 \times B$, and if that exceeds available memory the plan is rejected before any data is processed.

Numerical concept 4: the “grain” — how many elements per thread

Parallelizing has a fixed overhead: spawning a thread and coordinating it costs time. If a dataset is tiny, splitting it across threads costs more than it saves. So the resolver needs a grain: the minimum number of elements per thread that justifies a thread. The number of threads is ceil(elements / grain), clamped to the core count. A 100-element dataset below one grain (e.g. 1,024 elements) resolves to one thread; a 10-million-element dataset resolves to ceil(10,000,000 / 1,024) = 9,766, clamped to, say, 8 cores. This is why parallelism is part of the plan.

Observability concept: tracing spans

tracing is a Rust library for structured, hierarchical logging. A span is a named region of execution with a start, an end, and structured fields (key–value pairs). Spans can nest (a fit span containing a predict span). A span is created, its fields are filled in over time (some like duration_ms and error start empty and are filled at the end), and it is finalized when the operation completes. An OpenTelemetry layer is a piece that plugs a subscriber (the thing collecting spans) together with an exporter that sends them to a collector (like Jaeger). Making that layer opt-in behind a feature flag keeps the default build free of any network exporter.

4. Each object, explained

SKExecutionMode

What it is: A #[non_exhaustive] enum with five variants — Automatic, InProcessSynchronous, InProcessAsynchronous, OutOfCoreStreaming, OutOfCoreMemoryMapped. It is the intent the consumer declares; it is stored in the builder and can be Copy + Clone + PartialEq (cheap to duplicate and compare).

Why it exists: It separates “what the user wants” from “what will actually happen”. The user picks one value (or leaves it Automatic); the resolver later turns that into a plan. The Automatic variant is the key: it is the default that lets the library decide.

What we rejected: Free-form strings and allowing arbitrary extra modes without a fixed list. The enum makes invalid modes a compile-time error, and non_exhaustive still lets us grow the set later.

SKAccessPattern

What it is: A #[non_exhaustive] enum with three variants — Sequential, RandomAccess, Iterative. It is not a user choice; it is a property the algorithm declares about how it must touch its data.

Why it exists: It is the half of the compatibility check the user never controls. The resolver consults it to decide whether a requested mode is even possible: random access can never stream.

What we rejected: Letting the user set the access pattern (they would not know the algorithm’s internals) and deriving it heuristically at runtime (unsafe, opaque). The algorithm author declares it once, honestly, and the type keeps it honest.

SKExecutionContext

What it is: A plain struct holding the inputs to resolution: available_memory_bytes, cpu_cores, dataset_size_bytes, dataset_elements, scalar_size_bytes, grain, access_pattern, and batch_size_hint. It has two constructors: real() reads the physical machine (memory via sysinfo, cores via thread::available_parallelism), and tests build their own by hand with simulated values.

Why it exists: It makes resolution pure — all the environment enters through one explicit argument rather than hidden global state. That is what lets tests inject a fake 1-byte machine and get deterministic answers, and it is what makes “same inputs, same plan” a guarantee.

What we rejected: Reading the machine inside the resolver. That couples every resolution to the physical environment, breaks reproducibility, and makes the resolver impossible to unit test.

SKExecutionPlan

What it is: A struct with four fields — mode, parallelism, batch_size: Option<usize>, and buffer_depth. It is the output of resolution, PartialEq + Eq so two plans can be compared for equality.

Why it exists: It is the concrete, executable decision: “run in this mode, on this many threads, in batches of this size (if streaming), prefetching this many batches ahead.” It is what an algorithm actually consumes to run an operation.

What we rejected: Confusing the plan with the intent (they are deliberately different types), and leaving parallelism/batch_size out — an executable plan needs all of them.

sk_resolve_execution_plan

What it is: The free-scope public function sk_resolve_execution_plan(intent, context) -> Result<SKExecutionPlan, SKError>. It matches the intent; for Automatic it compares dataset_size_bytes to available_memory_bytes. If it fits, the mode is InProcessSynchronous; if not, the mode follows the access pattern (RandomAccess → OutOfCoreMemoryMapped; Sequential/Iterative → OutOfCoreStreaming). Explicit modes are kept unless a streaming intent meets a random-access pattern, which returns ExecutionModeIncompatible. Parallelism comes from parallelism_for(mode, context), and the buffer guard check_buffer_guard(context) rejects an oversized double-buffer footprint.

Why it exists: It is the pure, deterministic heart of the whole design (spec execution-decision). “Intent + context → plan” with no hidden state, so it is safe to call from any number of threads concurrently.

What we rejected: Resolving against the physical machine directly, and allowing automatic to ever produce an incompatibility error (it never does — automatic always picks a compatible mode).

sk_run_streaming_driver

What it is: The free-scope function that actually runs a streaming operation: sk_run_streaming_driver(source, state, update, plan) -> Result<(), SKError>. It spawns a dedicated producer thread that pulls batches from an SKLazySource into a bounded double buffer, while the caller’s thread consumes each owned SKDataBatch and applies the algorithm-supplied update(batch, parallelism, state) closure. update returns an SKStreamDecision (Continue or Stop) to steer the pipeline.

Why it exists: It implements the prefetch contract — input and computation overlap, batches arrive in order and fully owned, and a mid-stream read error surfaces as a central SKError with the effects of earlier batches preserved. This is the executor the plan’s streaming modes drive.

What we rejected: Reading and computing on the same thread (no overlap), and letting the algorithm push unboundedly ahead (hence the bounded buffer of buffer_depth).

SKBuilderState

What it is: A struct holding one field, execution_intent: SKExecutionMode, defaulting to Automatic. It offers execution_mode(...) (the chainable setter), execution_intent() (the getter), and resolve_plan(context) which delegates to sk_resolve_execution_plan.

Why it exists: It is the reusable “storage drawer” every algorithm builder embeds, so the execution-intent plumbing is written once instead of per algorithm (spec base-builders). It carries the Automatic default and knows how to turn stored intent into a concrete plan.

What we rejected: Each algorithm re-implementing intent storage and resolution — duplicated logic and drifting defaults. A single shared state type keeps the default uniform everywhere.

SKBuilder

What it is: A trait with two methods: execution_mode<E>(&mut self, mode: E) -> &mut Self where E: TryInto<SKExecutionMode> and build(self) -> Result<Model, SKError>.

Why it exists: It is the mandatory contract every estimator builder implements, so all builders look the same: configure execution, call build(), get a Result. It encodes “construction always goes through build()” at the type level (spec base-builders).

What we rejected: Letting algorithms skip the builder (they can’t — the trait and shared state are the foundation every one compiles against) and build() returning a bare model instead of a Result (loses the guarantee that invalid config surfaces as an error, not a panic).

SKOperationKind and SKBackendKind

What it is: Two enums of bounded strings. SKOperationKind is the closed set of library operations that emit spans (Fit, PartialFit, FitTransform, Transform, Predict, FitPredict, Score). SKBackendKind is the closed set of math backends (Faer, MatrixMultiply, NdArrayLinalg). Both implement Display so they render as clean strings like "fit" and "faer".

Why it exists: Structured observability needs bounded, compiler-checked values. Modelling these as enums means a misspelled operation or backend is a compile-time error, not a runtime surprise (spec observability: “the compiler rejects invalid values … instead of letting them surface at runtime”).

What we rejected: Free-form String attributes in the span — a typo would silently produce a mislabeled trace.

SKOperationAttributes

What it is: A struct bundling the structured fields of an operation span: operation (SKOperationKind), rows, columns, execution_mode (SKExecutionMode), and backend (SKBackendKind).

Why it exists: It is the single, type-safe bundle handed to the observability system when a span opens — one argument instead of five, all compiler-checked. Note it carries the resolved execution_mode, so a trace tells you which mode an operation actually ran in.

What we rejected: Passing each field separately to the span opener, and letting the span record raw, untyped values.

SKOperationObservation

What it is: A live span wrapper. It is created with begin(attributes), runs work via in_scope(...) (or enter()), records a failure with record_error(&error), and closes with finish(), which stamps the elapsed duration_ms onto the span.

Why it exists: It encapsulates the tracing::Span lifecycle so callers never touch tracing directly and never forget to fill duration_ms or record an error.

What we rejected: Exposing the raw tracing::Span and making callers manage its lifecycle — error-prone, and easy to leak a span.

sk_run_operation

What it is: A convenience function: sk_run_operation(attributes, operation) -> Result<T, E>. It opens an observation, runs the closure operation, records any error on the span, finishes the span, and returns the result unchanged.

Why it exists: It is the one-liner algorithm crates wrap each public operation in, so every fit/predict/transform is automatically traceable with failure visible in the trace (spec observability: “each public operation … so failures stay visible in the trace”).

What we rejected: Requiring each algorithm to hand-roll span open/error/close — duplication and easy to get wrong.

sk_opentelemetry_layer

What it is: A free-scope function that builds a Box<dyn tracing_subscriber::Layer<...>> wiring tracing to an OpenTelemetry tracer, given a tracer_name. It is only compiled when the observability-export feature is on.

Why it exists: It is the opt-in export path. The consumer configures the OpenTelemetry tracer provider and combines the returned layer with their subscriber; the default build compiles this out entirely and never initializes an exporter (spec observability: “the default build emits no OTLP traffic”).

What we rejected: Initializing an exporter by default, or coupling the library to a specific tracing subscriber — both would force telemetry behavior on consumers who don’t want it.

5. A real walk-through, using the tests

The most convincing way to learn is to read actual passing tests. Here are three — one from execution_tests.rs and two from streaming_driver_tests.rs.

Test 1: automatic_oversized_sequential_streams

This test proves that when an automatic intent meets a dataset larger than memory with a sequential access pattern, resolution picks streaming with the declared batch size.

#![allow(unused)]
fn main() {
#[test]
fn automatic_oversized_sequential_streams() {
    let ctx = simulated(
        1 << 30,          // available memory = 1 GiB
        4,                // CPU cores
        1 << 40,          // dataset size = 1 TiB  → larger than memory!
        None,             // dataset_elements unknown
        8,                // scalar size = 8 bytes (f64)
        SKAccessPattern::Sequential,
        Some(1024),       // batch size hint = 1024 elements
    );
    let plan = sk_resolve_execution_plan(SKExecutionMode::Automatic, &ctx).unwrap();
    assert_eq!(plan.mode, SKExecutionMode::OutOfCoreStreaming);
    assert_eq!(plan.batch_size, Some(1024));
}
}

Line by line:

1. let ctx = simulated(...) — we call the test helper simulated, which builds an SKExecutionContext by hand, never touching the physical machine. The 1 << 30 is bit-shift notation for 2^30 = 1 GiB of available memory; 1 << 40 is 2^40 = 1 TiB of dataset. The dataset is a thousand times bigger than memory.
2. SKAccessPattern::Sequential tells the resolver the algorithm reads in order, so streaming is a legal option. Some(1024) is the caller’s batch-size hint.
3. sk_resolve_execution_plan(SKExecutionMode::Automatic, &ctx) — automatic intent. Inside, the resolver sees dataset_size_bytes (1 TiB) <= available_memory_bytes (1 GiB) is false, so the dataset does not fit in memory. It then matches the access pattern Sequential and picks OutOfCoreStreaming. .unwrap() asserts the Result is Ok.
4. assert_eq!(plan.mode, SKExecutionMode::OutOfCoreStreaming) — confirms the mode.
5. assert_eq!(plan.batch_size, Some(1024)) — confirms the batch hint flowed into the plan, because in streaming mode batch_size is the hint (in other modes it is None).

The lesson: a single automatic intent produced a streaming plan purely from context — no user flag, no hidden machine reads.

Test 2: ordered_delivery_of_owned_batches

This test proves the streaming driver hands batches to the compute step in source order, and that each batch’s data survives even after the source has moved past it (i.e. the batch is owned, a copy, not a borrow).

#![allow(unused)]
fn main() {
#[test]
fn ordered_delivery_of_owned_batches() {
    let (mut source, _read_log, _outstanding, _max) = RecordingSource::new(5, None);
    let update =
        |batch: SKDataBatch<f64>, _parallelism: usize, state: &mut Vec<(usize, Vec<f64>)>| {
            let data = batch.data().as_slice_memory_order().unwrap().to_vec();
            state.push((batch.position(), data));
            SKStreamDecision::Continue
        };
    let mut state = Vec::new();
    sk_run_streaming_driver(&mut source, &mut state, update, &streaming_plan()).unwrap();

    assert_eq!(state.len(), 5);
    for (index, (position, data)) in state.iter().enumerate() {
        assert_eq!(*position, index);
        assert_eq!(data, &[index as f64]);
    }
}
}

Line by line:

1. RecordingSource::new(5, None) — the test’s instrumented source, total = 5 batches, fail_at = None (never fails). It also returns three shared log handles (read_log, outstanding, max) we ignore here.
2. The update closure is the algorithm’s per-batch step. Its signature is (batch, parallelism, state). Here state is a Vec<(position, data)> we accumulate. Inside, batch.data() returns a view of the batch’s matrix; .as_slice_memory_order() .unwrap().to_vec() copies it into an owned Vec<f64>, and we push (position, data). We return SKStreamDecision::Continue to keep going.
3. let mut state = Vec::new() — the algorithm’s mutable accumulator.
4. sk_run_streaming_driver(&mut source, &mut state, update, &streaming_plan()).unwrap() — run the pipeline. The driver’s producer thread pulls batches from source into the double buffer; the caller’s thread calls update for each in order. .unwrap() asserts success.
5. The assertions: state.len() == 5 (all five batches delivered), and for each index, the position equals the index and the data equals [index as f64] — i.e. batch k arrived at step k, holding the value k. Because we copied into an owned Vec, the assertion works even though the source has already produced the next batch — proving ownership.

The lesson: ordering is exact, and each batch is a fully owned SKDataBatch, safe to use after the source advances — exactly what the double-buffer prefetch promises.

Test 3: early_stop_halts_the_pipeline

This test proves an update returning Stop halts the pipeline, keeps the current batch’s effect, reads no batch beyond the one prefetched, and returns without error.

#![allow(unused)]
fn main() {
#[test]
fn early_stop_halts_the_pipeline() {
    let (mut source, read_log, _outstanding, _max) = RecordingSource::new(8, None);
    let update = |batch: SKDataBatch<f64>, _parallelism: usize, state: &mut Vec<usize>| {
        state.push(batch.position());
        if batch.position() == 2 {
            SKStreamDecision::Stop
        } else {
            SKStreamDecision::Continue
        }
    };
    let mut state = Vec::new();
    let result = sk_run_streaming_driver(&mut source, &mut state, update, &streaming_plan());

    assert!(result.is_ok(), "early stop should return without error");
    assert_eq!(state, vec![0, 1, 2]);
    let max_read = read_log
        .lock()
        .unwrap()
        .iter()
        .map(|(index, _)| *index)
        .max()
        .unwrap();
    assert!(
        max_read <= 3,
        "read beyond the prefetched batch: up to index {max_read}"
    );
}
}

Line by line:

1. RecordingSource::new(8, None) — a source of 8 batches, never failing. We keep read_log this time.
2. update pushes each batch’s position, then returns Stop exactly when position == 2, otherwise Continue.
3. sk_run_streaming_driver(...) — note we do not .unwrap(); we capture result.
4. assert!(result.is_ok()) — stopping early is a normal, successful end; no error.
5. assert_eq!(state, vec![0, 1, 2]) — the effect of batches 0, 1, 2 is kept; batch 2 was the last processed.
6. We inspect read_log (the recorded read timestamps). Since buffer_depth = 1 (double buffering), when compute stopped at batch 2 the producer had already prefetched batch 3 — so max_read <= 3 proves we read at most one batch ahead and then halted. No batch 4+.

The lesson: the driver hands control flow back to the algorithm through the return value of update — stopping is cooperative, keeps work already done, and bounds how far ahead the reader can go by the buffer depth.

6. Look inside

Here is the end-to-end flow from intent to streaming execution:

flowchart LR
    accTitle: From builder to streaming execution and observability
    accDescr: A builder holds the SKBuilderState with an execution intent. The resolver turns the intent and execution context into an SKExecutionPlan, or rejects incompatible or oversized plans. Streaming execution prefetches batches while updating state, and observability records each operation as an optional OpenTelemetry span.
    subgraph configure[Builder]
        BS[SKBuilderState] -->|intent<br/>SKExecutionMode::Automatic| RES
    end

    subgraph resolve[Resolution — pure & deterministic]
        CTX[SKExecutionContext<br/>memory, cores, size, pattern] --> RES
        RES[sk_resolve_execution_plan] --> PLAN[SKExecutionPlan<br/>mode · parallelism · batch_size · buffer_depth]
        RES -. incompatibility .-> ERR[SKError::ExecutionModeIncompatible]
        RES -. buffer overflow .-> ERR2[SKError::BatchBufferOversize]
    end

    PLAN -->|OutOfCoreStreaming| DRIVER
    subgraph run[Execution]
        DRIVER[sk_run_streaming_driver]
        SRC[SKLazySource] -->|batch k+1 prefetch| DRIVER
        DRIVER -->|owned SKDataBatch k| UPDATE[update batch → state]
        UPDATE -->|Continue / Stop| DRIVER
    end

    UPDATE --> OBS
    subgraph observe[Observability]
        OBS[sk_run_operation / SKOperationObservation<br/>span: operation · shape · mode · backend · duration · error]
        OBS -->|opt-in feature| OTL[sk_opentelemetry_layer → exporter]
    end

    style configure fill:#e0f2fe,stroke:#0284c7,color:#0c4a6e
    style resolve fill:#fef3c7,stroke:#d97706,color:#78350f
    style run fill:#f0fdf4,stroke:#16a34a,color:#14532d
    style observe fill:#ede9fe,stroke:#7c3aed,color:#3b0764
    style ERR fill:#ffe4e6,stroke:#e11d48,color:#7f1d1d
    style ERR2 fill:#ffe4e6,stroke:#e11d48,color:#7f1d1d

Reading the diagram: the builder holds the intent; resolution combines it with an injected context to produce a concrete plan (or a precise error); a streaming plan drives the double-buffered driver, which feeds owned batches to the algorithm’s update; and every operation runs inside an observability span, exported only if the feature is on.

7. Recap

• Intent and plan are different. The user stores a loose SKExecutionMode (defaulting to Automatic); the operation resolves it against a context into a concrete SKExecutionPlan with mode, parallelism, batch size, and buffer depth.
• Resolution is pure and deterministic. sk_resolve_execution_plan(intent, context) reads nothing global, so the same inputs always yield the same plan and tests can inject a fake machine.
• Automatic never fails; explicit errors are loud. Automatic always finds a compatible mode; an explicit streaming intent on a random-access algorithm fails before processing any data with ExecutionModeIncompatible, and an oversized double-buffer footprint fails with BatchBufferOversize.
• The builder pattern is mandatory. Every estimator goes through a typed builder embedding SKBuilderState and implementing SKBuilder, so construction always passes through build() → Result and the execution intent always defaults to Automatic.
• Streaming double-buffers. sk_run_streaming_driver overlaps input and compute through a bounded buffer, delivers ordered owned batches, and lets update stop the pipeline early.
• Observability is opt-in. Every operation runs under a tracing span with structured fields; exporting to OpenTelemetry is gated behind the observability-export feature so the default build emits no telemetry traffic.

Next chapter: now that an operation can resolve how it runs and be traced while it runs, we turn to the computational core — how the heavy math kernels are selected and parallelized under the hood.

The host-centric backend

The backend that does the heavy math started out married to one linear-algebra library (faer) and one numeric type (f64), so every algorithm that wanted a matrix product or a decomposition was stuck with that library’s data types. This chapter is the story of how we broke that marriage: we made the backend host-centric — it speaks only in the library’s own vocabulary (ndarray views in, ndarray results out) — and made the choice of who actually computes replaceable, without touching a single caller.

1. The problem we were solving

Imagine you run a restaurant. The kitchen uses one specific brand of stove, and every menu item is written as “cook it on that exact stove.” Now you want to add a faster stove, or a cheaper one, and offer both. You would have to rewrite every recipe. Worse, the waiters (your algorithms) would have to learn how each stove works just to write down an order.

That was our situation. The math backend — the code that multiplies matrices and factorizes them — used faer, a pure-Rust linear-algebra library, and hard-coded the floating-point type f64. Its public methods took faer::Mat and returned faer::Mat. So every algorithm that wanted a matrix operation had to think in faer’s types. There was no way to swap in a system-optimized BLAS library (OpenBLAS, MKL, Accelerate) without rewriting the algorithms. And there was only one numeric type, f64 — no room for the lighter, faster f32 that many models use to save memory.

What goes wrong without a fix? Every new algorithm (Linear Regression, PCA, Ridge, GMM — the roadmap was filling up) would be welded to faer, and adding a second backend would mean touching everything. The library could not grow. So we asked a simple question: what is the smallest, cleanest promise a math backend can make to the rest of the library?

The answer became the guiding idea of this whole chapter: the backend adapts to the library, not the other way around. The library already has its own container type, ndarray, for dense numbers. So the backend’s public surface is written in ndarray terms — you hand it a view of your numbers, it hands back an owned result. Whether faer, a system BLAS, or (much later) a GPU sits behind the curtain is the backend’s private business.

2. The decisions — and the roads not taken

3. The concepts, taught from zero

Before we read any code, let’s build the two ideas you need: what a trait is, and what “host-centric” really means. We’ll also touch on the small vocabulary the code uses.

A trait is a contract

A trait in Rust is a set of promises a type makes. Think of an electrical outlet: the outlet promises “any device with the right plug can draw power here.” The device doesn’t care where the electricity comes from, and the outlet doesn’t care what the device is. A trait works the same way.

#![allow(unused)]
fn main() {
pub trait SKMathBackend<F: SKFloat>: Send + Sync {
    fn kind(&self) -> SKBackendKind;
    fn gemm(&self, a: ArrayView2<F>, b: ArrayView2<F>, alpha: F, parallelism: usize) -> Array2<F>;
    // ... and many more
}
}

This says: “any type that wants to be a math backend must provide a kind, a gemm, a svd, and so on.” The outlet (the trait) doesn’t say how to compute — it just says the plug fits. A backend like SKFaerBackend then declares “I am a math backend” and fills in each method using its own library.

Generic over <F: SKFloat>

The <F: SKFloat> is a generic parameter — a placeholder for a type, the way a recipe’s “add flour” leaves how much flour open. The little : SKFloat is a bound: it means F is not any type, only types that implement SKFloat. And SKFloat is sealed: only f32 and f64 can ever be SKFloat, because only those two are declared to implement it. So <F: SKFloat> is a polite way of saying “this works for f32 and f64, and nothing else.”

The whole trait is written once, generically over F. Then for each concrete float we plug in the real types — ArrayView2<f64>, ArrayView2<f32>, and so on — automatically.

dyn vs generics: why a box, not a concrete type

Rust has two ways to use a trait. Generics bake the concrete type in at compile time: SKFaerBackend::new().gemm(...) knows exactly which type it is, so the compiler can make everything fast and type-safe. But generics make the caller name the type.

Sometimes the caller doesn’t know — or doesn’t care — which backend it has. It just wants “a math backend.” That’s where dyn comes in. Box<dyn SKMathBackend<F>> is a trait object: a pointer to some value that happens to implement SKMathBackend<F>. The caller holds the box and calls the trait’s methods; the actual type is decided at runtime (the dispatch function decides, below).

Why not always use dyn? Because a trait object adds a small runtime indirection and the compiler can’t inline it. Why not always use generics? Because the default backend is chosen by build configuration, not by the caller — we literally do not know which type will be there. So we use generics for the concrete backends (fast, in algorithms that name them) and dyn at the boundary where configuration picks the type (Box<dyn SKMathBackend<F>>). Each tool for its job.

What “host-centric” means

The word host means “the machine the program runs on” — your CPU and RAM, as opposed to a separate device like a GPU. Host-centric is a design rule: the library’s public vocabulary is the host’s vocabulary (ndarray for dense data, sprs for sparse data), and the backend must translate to and from it at the door.

Concretely: the trait’s inputs are ArrayView2<F> — views, which are a borrowed window over data that already exists (zero copies). Its outputs are Array2<F> — owned containers, fresh memory the caller takes possession of. No faer type appears anywhere in the public surface. If we swap the backend from faer to ndarray-linalg, every algorithm keeps talking to ndarray exactly as before; only the door changes.

#![allow(unused)]
fn main() {
fn assert_host_centric_surface<F: SKFloat, B: SKMathBackend<F> + ?Sized>(
    backend: &B,
    a: Array2<F>,
) {
    let _ = backend.gemm(a.view(), a.view(), F::one(), 1);
    let _: Result<SKSingularValueDecomposition<F>, _> = backend.svd(a.view());
    let _: SKQRDecomposition<F> = backend.qr(a.view());
    // ...
}
}

This test compiles if and only if the trait’s public surface is free of faer types. It’s the compiler enforcing “host-centric” for us — a theme we’ll revisit in the walk-through.

A small vocabulary of the code

• ndarray — the library’s dense-array container. Array2<F> is a 2-D grid of F (a matrix); ArrayView2<F> is a read-only view over one (borrowed, no copy).
• SKFloat — the sealed scalar bound; only f32 and f64.
• SKBackendKind — a small enum naming which backend produced a result (Faer, MatrixMultiply, NdArrayLinalg), recorded in observability spans.
• SKError — the crate’s central error type; backends map their library’s errors into it so callers never see a faer or LAPACK error.
• Decompositions — ways to rewrite a matrix into simpler pieces: SVD, QR, LU, Cholesky. They’re the building blocks of many ML algorithms.

4. Each object, explained

SKMathBackend<F> (trait)

What it is — The single contract for “a thing that can do heavy dense algebra.” It declares the operations the library needs: gemm (matrix multiply), svd, qr, cholesky, solve_triangular, solve, eigh, lu, slogdet, pinv, inv, lstsq, norm, and vector_norm. It is generic over F: SKFloat, so the same contract serves both f32 and f64. It is Send + Sync, so a math operation can be shared across threads (an execution plan can fan the work out).

Why it exists — It is the abstraction boundary. Everything above it (the algorithms, the execution planner) talks only to this trait, in ndarray terms. Everything below it (the concrete backends) does the actual number crunching. The library is written once, against the trait, and any backend that honors the trait plugs in.

What we rejected — We rejected leaking concrete types into the trait (the old f64 + faer::Mat surface), associated matrix types that would make the surface device-ready now, and making the trait async (the backend stays eager and pure; concurrency is the execution planner’s job). Also rejected: hard-coding a sequential GEMM — the trait takes a parallelism argument so large products scale with the plan.

SKNormKind<F> (enum)

What it is — The set of matrix/vector norm orders the backend understands: Frobenius, L2, L1, Infinity, their negatives, Nuclear, and a General { order: F } fallback for arbitrary p-norms. It mirrors numpy.linalg.norm / scipy.linalg.norm.

Why it exists — Many algorithms need a norm (a single number measuring “size”) in different flavors. Instead of a zoo of separate methods, one norm / vector_norm method takes an SKNormKind to say which flavor.

What we rejected — We rejected a separate method per order, and we kept the General p-norm fallback rather than only hard-coding the common cases.

The decomposition containers

These four structs are the results a backend returns. Crucially, each holds plain Array2<F> fields — owned ndarray matrices — so the result does not depend on which backend produced it. That is the host-centric promise made visible: you can take a decomposition from faer, hand it to the BLAS backend, and both are the same type.

SKSingularValueDecomposition<F>

What it is — The result of a singular value decomposition $A = U \Sigma V^{\mathsf{T}}$: fields u (left singular vectors, Array2<F>), singular_values (a Vec<F>, non-increasing), and v (right singular vectors, Array2<F>).

Why it exists — SVD is the workhorse behind pseudo-inverses, lstsq, PCA, and the spectral norm. A concrete, backend-independent container means an algorithm can use the result without caring who computed it.

What we rejected — Backend-native SVD types (faer’s Svd, LAPACK’s), which would leak the backend into every consumer.

SKQRDecomposition<F>

What it is — The result of a (thin) QR decomposition A = Q R: fields q (orthogonal factor) and r (upper-triangular factor), both Array2<F>.

Why it exists — QR is a numerically stable way to solve least-squares problems and to find orthonormal bases.

What we rejected — Any backend-specific QR representation.

SKLUDecomposition<F>

What it is — The result of an LU decomposition P A = L U: fields pivot (a host Vec<usize> row permutation where pivot[i] is the original row that lands at position i), l (unit-lower-triangular), and u (upper-triangular).

Why it exists — LU is the engine behind solve and inv. The pivot is kept as a plain Vec<usize> so the row permutation is portable and readable.

What we rejected — A PivotInfo associated type (design Decision 2): every result is host-resident, so there is nothing backend-specific to abstract.

SKLeastSquaresSolution<F>

What it is — The result of a minimum-norm least-squares solve $x = \operatorname{argmin} , |b - A x|_2$: fields solution (Array2<F>), rank (the effective rank after rcond truncation), singular_values (Vec<F>), and residual_sum_of_squares (an Option<Vec<F>>, present when the system is overdetermined, i.e. more rows than columns).

Why it exists — It bundles everything an algorithm (like Linear Regression) needs to report how well the fit worked, in one backend-neutral struct.

What we rejected — Returning only the solution and making callers recompute rank and residuals — duplicated tolerance logic per caller.

SKFaerBackend

What it is — The pure-Rust default backend. A tiny, zero-field struct that implements SKMathBackend<F> for both f32 and f64 by delegating to the faer library. GEMM uses faer’s blocked matmul (honoring parallelism); decompositions delegate to faer’s high-level SVD/QR/Cholesky/LU/eigendecomposition. No C or Fortran BLAS is linked.

Why it exists — It is the out-of-the-box path: compiles on Rust 1.85 with no external toolchain, so a default build “just works” for everyone (spec blas-interface).

What we rejected — Making a C BLAS the default (would force toolchain setup); using oxiblas (rejected because it fails to compile on Rust 1.85).

SKMatrixMultiplyBackend

What it is — A second pure-Rust backend that routes heavy GEMM to the matrixmultiply crate’s dense kernel. Decompositions aren’t in matrixmultiply, so they delegate to faer, which is always available.

Why it exists — It gives an explicit, non-faer GEMM path — useful as an experiment or fallback for matrix products specifically.

What we rejected — Nothing major; it’s a lean companion to the faer backend, and it exercises the same SKMathBackend contract.

SKNdArrayLinalgBackend

What it is — The opt-in system-BLAS backend, backed by ndarray-linalg over a configured blas-src library (OpenBLAS, BLIS, MKL, or Accelerate). It implements the same SKMathBackend<F> surface for f32 and f64. Its lstsq uses the native LAPACK gelsd driver, and its solve iterates columns of b using the native solve.

Why it exists — When you want the peak performance of a hand-tuned BLAS on your machine, you enable the blas-backend feature and this backend takes over heavy dense algebra. It is gated behind the blas-backend feature and is never on by default.

What we rejected — Making it the default (toolchain cost), and duplicating the whole feature behind an unconditional dependency.

sk_default_math_backend<F>()

What it is — A free function that returns the backend configured for the current build: Box<dyn SKMathBackend<F>>. With blas-backend enabled it returns SKNdArrayLinalgBackend; otherwise SKFaerBackend.

#![allow(unused)]
fn main() {
#[cfg(not(feature = "blas-backend"))]
pub fn sk_default_math_backend<F: SKFloat>() -> Box<dyn SKMathBackend<F>>
where
    SKFaerBackend: SKMathBackend<F>,
{
    Box::new(SKFaerBackend::new())
}
}

Why it exists — The execution planner needs “a math backend” but should not care which build was chosen. This is the single place where build configuration becomes a concrete backend. The #[cfg(...)] attributes mean the compiler picks the right body based on the blas-backend feature — a compile-time decision, so there is no runtime cost to choosing.

What we rejected — Baking the choice into every call site, or forcing callers to name a concrete type (they can’t — they don’t know the build config).

The input wrapper: FaerInput (private)

What it is — A small private enum inside faer_backend.rs that wraps a faer matrix view: either Borrowed(MatRef<'a, F>) for zero-copy input, or Owned(Mat<F>) when a copy was needed.

Why it exists — This is where “host-centric” meets “zero-copy” in practice. When the ndarray input is C-contiguous (row-major), the backend builds a faer view over the same memory — no copy. When the input has a non-standard (strided) layout, faer needs contiguous data, so the backend makes an internal copy. The public surface is always a zero-copy view; the copy, when it happens, is the backend’s private concern (design’s risk note).

What we rejected — Forcing every caller to guarantee contiguous memory (too strict), or always copying (wastes the zero-copy benefit). The Borrowed/Owned pair gives the best of both.

#![allow(unused)]
fn main() {
enum FaerInput<'a, F: SKFloat> {
    Borrowed(MatRef<'a, F>),
    Owned(Mat<F>),
}
}

5. A real walk-through, using the tests

Let’s read three tests from backend_tests.rs and see the ideas in action. The whole suite is deliberately host-centric: fixtures are built from ndarray, and no faer type appears in any test — that’s the discipline the design enforces.

Test 1: GEMM matches a hand-computed reference product

#![allow(unused)]
fn main() {
fn check_gemm_matches_reference_product<F: SKFloat + std::fmt::Debug>()
where
    SKFaerBackend: SKMathBackend<F>,
{
    let a = mat_from_rows(&[
        &[F::one(), F::from(2.0).unwrap()],
        &[F::from(3.0).unwrap(), F::from(4.0).unwrap()],
    ]);
    let b = mat_from_rows(&[
        &[F::from(5.0).unwrap(), F::from(6.0).unwrap()],
        &[F::from(7.0).unwrap(), F::from(8.0).unwrap()],
    ]);
    let backend = SKFaerBackend::new();
    let product = backend.gemm(a.view(), b.view(), F::one(), 1);
    let expected = mat_from_rows(&[
        &[F::from(19.0).unwrap(), F::from(22.0).unwrap()],
        &[F::from(43.0).unwrap(), F::from(50.0).unwrap()],
    ]);
    assert_close(&product, &expected, F::from(1e-4).unwrap());
}
}

Line by line:

• The function is generic over F and demands SKFaerBackend: SKMathBackend<F> — i.e. “run this for whatever float makes the backend valid.” The same check then runs for f64 and f32 via two tiny #[test] wrappers.
• mat_from_rows builds a 2-D Array2<F> from slices. F::one() and F::from(2.0) make values without assuming a concrete float type — this is how you write “the number two” generically.
• We multiply a (the matrix [[1,2],[3,4]]) by b ([[5,6],[7,8]]) with alpha = 1.0. The 1 is parallelism (sequential).
• The expected matrix is computed by hand: $1 \cdot 5 + 2 \cdot 7 = 19$, $1 \cdot 6 + 2 \cdot 8 = 22$, $3 \cdot 5 + 4 \cdot 7 = 43$, $3 \cdot 6 + 4 \cdot 8 = 50$. This is the “reference” — what matrix multiply must produce, no matter which backend.
• assert_close checks every element is within 1e-4 (floating-point math is never exact, so we allow a tiny tolerance). The test proves the backend’s gemm is correct, not just that it runs.

Test 2: the trait surface stays host-centric (the compiler is the referee)

#![allow(unused)]
fn main() {
fn assert_host_centric_surface<F: SKFloat, B: SKMathBackend<F> + ?Sized>(
    backend: &B,
    a: Array2<F>,
) {
    let _ = backend.gemm(a.view(), a.view(), F::one(), 1);
    let _: Result<SKSingularValueDecomposition<F>, _> = backend.svd(a.view());
    let _: SKQRDecomposition<F> = backend.qr(a.view());
    let _ = backend.cholesky(a.view());
    let _ = backend.solve_triangular(a.view(), a.view(), true, false);
    let _ = backend.solve(a.view(), a.view());
    let _ = backend.eigh(a.view());
    let _ = backend.lu(a.view());
    let _ = backend.slogdet(a.view());
    let _ = backend.pinv(a.view());
    let _ = backend.inv(a.view());
    let _ = backend.lstsq(a.view(), a.view());
    let _ = backend.norm(a.view(), SKNormKind::Frobenius);
    let _ = backend.vector_norm(a.slice(ndarray::s![.., 0]), SKNormKind::L2);
}
}

Line by line:

• This function is generic over any backend B: SKMathBackend<F> (note ?Sized, which allows dyn types too). It doesn’t build a backend or check results — it only exercises the whole public surface with ndarray views.
• Every call passes ArrayView2<F> (a.view()) and every result is typed as an Array2<F>-based container. Notice how we annotate the types (: Result<SKSingularValueDecomposition<F>, _>) — these annotations are the assertion. If the trait leaked a faer type anywhere, the function would not compile with a generic B.
• The test trait_surface_is_host_centric calls it with SKFaerBackend and with the default backend box:

  #![allow(unused)]
  fn main() {
  let backend = sk_default_math_backend::<f64>();
  assert_host_centric_surface(&*backend, a);
  }

  Both must compile — proving a dyn backend and a concrete backend present the same host-centric face. This is the design Decision 1 made executable: the compiler itself guarantees no backend type leaks.

Test 3: solve_triangular recovers a known solution

#![allow(unused)]
fn main() {
fn check_solve_triangular_upper<F: SKFloat + std::fmt::Debug>()
where
    SKFaerBackend: SKMathBackend<F>,
{
    let a = mat_from_rows(&[
        &[F::from(2.0).unwrap(), F::from(1.0).unwrap()],
        &[F::zero(), F::from(3.0).unwrap()],
    ]);
    let b = mat_from_rows(&[&[F::from(5.0).unwrap()], &[F::from(6.0).unwrap()]]);
    let backend = SKFaerBackend::new();
    let x = backend
        .solve_triangular(a.view(), b.view(), false, false)
        .unwrap();
    // 2x + y = 5; 3y = 6 => y = 2, x = 1.5
    assert!((x[(0, 0)] - F::from(1.5).unwrap()).abs() < F::from(1e-4).unwrap());
    assert!((x[(1, 0)] - F::from(2.0).unwrap()).abs() < F::from(1e-4).unwrap());
}
}

Line by line:

• a is upper-triangular: [[2,1],[0,3]] (everything below the diagonal is zero). We pass (false, false) — lower = false (it’s upper), unit_diagonal = false.
• We’re solving $A x = b$, i.e. $2x + y = 5$ and $3y = 6$. By hand: from the second equation, $y = 2$; plugging in, $2x + 2 = 5$, so $x = 1.5$. The comment spells this out.
• backend.solve_triangular(...).unwrap() asks the backend to solve it; .unwrap() unpacks the Result and panics if there was an error (fine for a test).
• Two assert!s check $x_{0,0} \approx 1.5$ and $x_{1,0} \approx 2.0$. Note we read the vector using two indices (row, col) — even a column vector is a 2-D Array2 here, so the single column is (0,0) and (1,0). This test checks numerical correctness of a less obvious operation, again generically over f64/f32.

Together, these three tests show the three pillars: correctness (test 1), a clean, backend-agnostic surface enforced by the compiler (test 2), and the expanded kernel working (test 3).

6. Look inside

Here is the host-centric boundary as a picture: the library speaks ndarray/sprs, and each backend adapts at the door, translating to and from its own library internally.

flowchart LR
    accTitle: The host-centric backend boundary
    accDescr: The library speaks host types such as ndarray views. Those cross the SKMathBackend trait, which takes array views in and returns owned arrays or decomposition containers out. The default SKFaerBackend adapts to faer, and the opt-in SKNdArrayLinalgBackend adapts to ndarray-linalg over a system BLAS.
    subgraph Host["The library speaks host types (ndarray / sprs)"]
        A["Algorithm (LinearRegression, PCA, ...)"]
        P["Execution planner"]
        C["SKDataView (ndarray views)"]
    end

    subgraph Boundary["SKMathBackend<F: SKFloat>  (the trait)"]
        I["ArrayView2 / ArrayView1 in"]
        O["Array2 / decomposition containers out"]
    end

    subgraph Faer["SKFaerBackend (default, pure Rust)"]
        F1["faer::Mat & MatRef"]
        F2["matrixmultiply GEMM"]
    end

    subgraph Blas["SKNdArrayLinalgBackend (opt-in blas-backend)"]
        N1["ndarray-linalg + blas-src"]
        N2["OpenBLAS / MKL / BLIS / Accelerate"]
    end

    A --> P
    C --> I
    P --> I
    I --> Boundary
    Boundary --> O
    O --> A
    O --> P

    Boundary -- "faer" --> Faer
    Boundary -- "system BLAS" --> Blas

    F1 -. "zero-copy view (contiguous)" .-> F2
    F2 --> O
    N1 --> O

The important detail: the arrows from the trait to Faer and Blas are private — the host never reaches past the boundary. The FaerInput wrapper is the door for faer: borrow the memory when the layout allows, copy only when it must.

7. Recap

• The backend is host-centric: it accepts ndarray views and returns owned ndarray results, so no concrete linear-algebra library leaks into the public surface.
• Generic over SKFloat, the sealed bound for f32/f64, gives every operation both float types with code written once.
• dyn vs generics is a judgment call: concrete backends are generic (fast), but the build-configured default is a Box<dyn SKMathBackend<F>> because the caller can’t know which type is active.
• faer is the pure-Rust default (works with no C toolchain); system BLAS (ndarray-linalg) is an opt-in feature, chosen at compile time by sk_default_math_backend.
• Decomposition containers are plain Array2<F> structs, so a result is portable across backends; zero-copy is preserved at the boundary, with internal copies only when input layout demands it.

Next, we’ll look at how the execution planner decides where and how to run work — the layer that calls this backend and composes the algorithms above it.

Module hygiene: folder modules & the tests-at-end rule

A machine-learning library starts as a handful of files and quietly becomes hundreds. Long before the algorithms get hard, the organization gets hard: you open a file to change one small thing and have to scroll through two thousand lines to find it. In this chapter we set up the discipline that keeps a growing Rust codebase navigable — a set of house rules about where code lives, not what it computes. We learn how a Rust file becomes a folder module, why a mod.rs should be a pure dispatcher that carries no logic, and why tests sit beside the implementation they check rather than being sprinkled through it or parked in a far-away tests/ directory.

1. The problem we were solving

Imagine a library that is really a filing cabinet. Early on, when the cabinet has three drawers, it does not matter much how you label them. But as the library grows — dozens of types, dozens of algorithms, thousands of lines — the cabinet fills up, and the way you organize it becomes the difference between finding a document in seconds or in hours.

sciencekit had reached that point. It had settled on a convention: when a source file grows past roughly 200 lines, you split it into a folder — a directory with a special mod.rs file that acts as the “label on the drawer”, telling the compiler what lives inside. And tests were supposed to live in companion files right beside the implementation they check.

But rules only help if every drawer follows them. Fourteen mod.rs files had drifted: instead of being just a label (a pure dispatcher), each of them had grown the actual implementation inline — the drawer label was also doing the filing. And there was no written rule at all about where inline tests must sit within a file. The result was a library where the same file had two jobs, where a small file became a big one and then a folder but nobody had finished the move, and where tests appeared at the top, middle, or bottom depending on whoever wrote them.

The change we walk through in this chapter, module hygiene, was a pure refactor: no new behavior, no new math, no changed results. It simply made every mod.rs a true pure dispatcher and codified the tests-at-end rule. It is the quiet maintenance work that keeps a codebase pleasant to live in, and it is a great lens for learning how Rust organizes code.

2. The decisions — and the roads not taken

Scope note: one of the fourteen modules, backend/, was deliberately excluded from this refactor. It was owned by a companion change (backend-kernel-expansion) that was restructuring it anyway. Splitting it twice would have caused conflicts, so this change covered the other thirteen and left backend/ alone.

3. The concepts, taught from zero

Rust: what a “module” is

In many languages you group related code into “files” and “packages”. Rust gives you a cleaner idea called a module. A module is simply a named container of code: it can hold functions, types, other modules, and even constants. You create one with the keyword mod:

#![allow(unused)]
fn main() {
mod shapes {
    pub struct Circle { radius: f64 }
}
}

This module is named shapes and contains a type Circle. The pub in front of struct says “let other code outside this module see it”. Without pub, the item stays private — visible only inside shapes.

Modules nest. A module inside a module, inside a module… This nesting is exactly what lets a big codebase be organized like a filing cabinet: each drawer is a module, and a drawer can contain smaller drawers.

How a file becomes a module, and then a folder module

Rust has a convenient rule: the file system can stand in for the module tree. If you write mod shapes; (note the semicolon, not a block), Rust looks for a file named shapes.rs and makes it the body of the shapes module.

But what if shapes grows too big for one file? Then you turn it into a folder:

shapes/            ← a directory
  mod.rs           ← the "label" file: declares the module's contents
  circle.rs        ← one implementation file
  rectangle.rs     ← another implementation file

The special file mod.rs inside a folder is the module. Rust reads mod shapes; and finds the shapes/ folder with its mod.rs, and everything inside that folder becomes part of the shapes module. This is the folder module pattern.

mod.rs as a pure dispatcher

Here is the key discipline of this whole chapter: mod.rs should contain no implementation logic at all. It has exactly two jobs:

1. Declare the submodules with mod <name>;.
2. Re-export their public items with pub use, so that code using the folder can reach the items without knowing which file they came from.

A mod.rs that only does these two things is called a pure dispatcher — it does no work itself; it just points to the files that do the work.

pub use: re-exports

pub use is a way to “re-share” an item. Suppose circle.rs defines Circle. Inside shapes/mod.rs you write pub use circle::Circle;. Now anyone who uses the shapes module can write shapes::Circle even though Circle physically lives in circle.rs. This gives you freedom: you can rearrange the files inside the folder without changing anything for the code that uses the folder. The mod.rs is the stable public face; the file layout is an internal detail.

#[cfg(test)]: the compile-time “only when testing” gate

#[cfg(...)] means “configuration”: it tells the compiler to include this item only under certain conditions. #[cfg(test)] means “only include this module when compiling in test mode”. So you can write:

#![allow(unused)]
fn main() {
#[cfg(test)]
mod shapes_tests;
}

This module is compiled only when tests run; in a normal release build it simply does not exist. This is how test code stays out of your shipped library.

Two ways to attach tests

You have two main choices for where test code lives:

• Inline: a #[cfg(test)] mod … block at the end of the same file as the implementation.
• Companion module: a separate *_tests.rs file (or folder) beside the implementation file, declared at the end of mod.rs.

sciencekit’s rule: tests live in companion *_tests.rs modules beside the implementation, and any inline tests must sit at the very end of their file. Never a global tests/ directory.

4. Each object, explained

The folder module (mod.rs + submodule files)

What it is: A directory whose mod.rs declares and re-exports its submodules, with each implementation living in its own file inside the folder. Triggered by the ~200-line rule: a single .rs file that crosses that length becomes a folder.

Why it exists: It bounds the size of any one file, gives every piece of code one obvious home, and keeps the public interface stable while the internals are free to rearrange. Opening mod.rs shows the whole catalogue at a glance.

What we rejected: Keeping long files as single giant .rs files (unreadable), and splitting into one file per function (over-fragmented). The folder module is the middle path: big enough to hold a real concern, small enough to read.

Let us look at two real pure-dispatcher mod.rs files. First, kernels from sciencekit_math — the folder that holds the higher-order math kernels:

#![allow(unused)]
fn main() {
//! Higher-order operation kernels built on `azip!`/`par_azip!`/`zip_mut_with`.

mod elementwise;
mod grain;
mod reductions;
mod scaling;

pub use elementwise::{sk_binary_combine, sk_elementwise_transform};
pub use grain::{
    SK_AXIS_SUM_GRAIN, SK_BINARY_COMBINE_GRAIN, SK_ELEMENTWISE_TRANSFORM_GRAIN,
    SK_SCALE_IN_PLACE_GRAIN,
};
pub use reductions::sk_axis_sum;
pub use scaling::sk_scale_in_place;

#[cfg(test)]
mod grain_tests;
#[cfg(test)]
mod kernels_tests;
}

See the structure? Four mod declarations, four pub use re-exports, then the #[cfg(test)] companion test modules at the end. Zero math logic here. If you want to know everything the kernels folder offers, you read this 27-line file and you are done.

Now execution from sciencekit_common — the folder that turns a user’s requested execution intent into a concrete plan:

#![allow(unused)]
fn main() {
//! Execution planning.

mod context;
mod modes;
mod plan;
mod resolve;
mod streaming_driver;

pub use context::SKExecutionContext;
pub use modes::{SKAccessPattern, SKExecutionMode};
pub use plan::SKExecutionPlan;
pub use resolve::sk_resolve_execution_plan;
pub use streaming_driver::{SKStreamDecision, sk_run_streaming_driver};

#[cfg(test)]
mod execution_tests;
#[cfg(test)]
mod streaming_driver_tests;
}

Same shape, slightly bigger: five submodules, five re-export lines, two companion test modules. The data_view and builders folders follow the identical pattern. This consistency is the point — once you have seen one pure-dispatcher mod.rs, you have seen them all.

The companion *_tests.rs module

What it is: A test module that lives beside the implementation it checks, declared with #[cfg(test)] mod …_tests; at the end of mod.rs. It may itself be a single file or a small folder of test files.

Why it exists: Tests and the code they verify are updated together, found together, and reasoned about together. Keeping them adjacent (but in a separate file) means the implementation file reads top-down while the tests stay close by, never in a distant global tests/ directory that drifts out of sync.

What we rejected: A global tests/ directory, and scattering #[cfg(test)] blocks throughout implementation files. Both separate tests from what they test or bury the implementation.

Here is a real companion test folder. The scorer_traits folder holds a mod.rs plus a scorer_traits_tests/ folder whose own mod.rs is:

#![allow(unused)]
fn main() {
//! Companion tests for the scoring contracts (split into a folder module to
//! respect the 200-line limit).
pub mod example_contracts;
pub mod scoring_tests;
}

It is itself a tiny pure dispatcher: it just declares the two test files inside it. This is the folder-module pattern applied recursively to the tests themselves — the tests grew past 200 lines, so they became a folder too.

Compare that with execution, where the two test modules are plain single files (execution_tests.rs, streaming_driver_tests.rs) because each stayed under the limit. The rule adapts: small tests, one file; big tests, a folder.

The 200-line rule

What it is: A stated convention that a .rs file over roughly 200 lines becomes a standardized folder module, with files named by role (mod.rs, core_implementation.rs, *_tests.rs, and for builders builder.rs).

Why it exists: It is a mechanical tripwire that triggers exactly when a file is doing too much. You do not need to argue about taste — you count lines, and past the limit you split. It keeps every file a comfortable size and gives every concern a named home.

What we rejected: No limit (files grow into impenetrable walls of text) and an emotionally enforced limit (inconsistent — one developer’s “fine” is another’s “too long”). A hard, stated number is uniformly applied.

5. A real walk-through, using the tests

Tests are our executable documentation, so let us read two real tests and see what the test-at-end, companion-module organization enables.

Test 1: prefetching under a streaming plan

This test lives in execution/streaming_driver_tests.rs — a companion file to the streaming_driver.rs implementation. It checks a subtle concurrency promise: while one batch is computing, the driver should already be reading the next one (this is called prefetch). The test is reader_advances_while_compute_runs:

#![allow(unused)]
fn main() {
#[test]
fn reader_advances_while_compute_runs() {
    let (mut source, read_log, _outstanding, _max) = RecordingSource::new(5, None);
    let finish_log: Arc<Mutex<Vec<(usize, Instant)>>> = Arc::new(Mutex::new(Vec::new()));
    let finish = Arc::clone(&finish_log);
    let update = move |batch: SKDataBatch<f64>, _parallelism: usize, _state: &mut ()| {
        thread::sleep(Duration::from_millis(50));
        finish
            .lock()
            .unwrap()
            .push((batch.position(), Instant::now()));
        SKStreamDecision::Continue
    };
    let mut state = ();
    sk_run_streaming_driver(&mut source, &mut state, update, &streaming_plan()).unwrap();

    let reads = read_log.lock().unwrap().clone();
    let finishes = finish_log.lock().unwrap().clone();
    let read_time = |index: usize| reads.iter().find(|(i, _)| *i == index).map(|(_, t)| *t);
    let finish_time = |index: usize| finishes.iter().find(|(i, _)| *i == index).map(|(_, t)| *t);

    // For every k except the last, batch k+1 was read before batch k finished.
    for k in 0..4 {
        let read_next = read_time(k + 1).expect("batch k+1 must be read");
        let done = finish_time(k).expect("batch k must complete");
        assert!(
            read_next < done,
            "batch {} read after batch {} finished — no prefetch overlap",
            k + 1,
            k
        );
    }
}
}

Let us read it slowly. The test creates a RecordingSource with 5 batches (the RecordingSource is a small instrumented stand-in defined right in the test file: it logs every read and records the time). The update closure is the “computation”: it simulates work by sleeping 50 milliseconds, then logs which batch finished and when. We then run the streaming driver over the source with this callback.

After the run, the test grabs two logs: reads (when each batch was read) and finishes (when each batch finished). The clever part is the loop. For every batch k from 0 to 3, it asserts that batch k + 1 was read before batch k finished:

#![allow(unused)]
fn main() {
assert!(read_next < done, "batch {} read after batch {} finished — no prefetch overlap", ...);
}

If the driver were naive and processed one batch fully before touching the next, batch k+1 would be read after batch k finished, and this assertion would fail. So this one test pins down the whole prefetch contract: the reader runs ahead while compute runs. That is a subtle piece of concurrency, expressed clearly in ~20 lines, sitting in a file right beside the streaming_driver.rs code it validates.

Test 2: pure scoring must not re-infer

This test lives deeper — inside the scorer_traits/scorer_traits_tests/scoring_tests.rs folder. It checks a design contract: a scorer has a “pure” form that scores from predictions you already have, and it must not run the model again. Here it is:

#![allow(unused)]
fn main() {
#[test]
fn label_pure_form_does_not_re_infer() {
    let scorer = accuracy::<f64, ConstantClassifier>();
    let truth = array![1_i64, 1, 0];
    let labels = array![1_i64, 1, 1];
    let score = <Accuracy as SKLabelScorer<f64, ConstantClassifier>>::score_from_labels(
        &scorer,
        SKTargetView::try_from(truth.view()).unwrap(),
        labels.view(),
    )
    .unwrap();
    // 2 of 3 correct (predicted 1,1,1 vs true 1,1,0).
    assert!((score - 2.0 / 3.0).abs() < 1e-9);
}
}

The test sets up truth [1,1,0] and labels [1,1,1]. A ConstantClassifier always predicts label 1, so 2 of the 3 labels match the true labels. The score should be 2/3. We call the pure form directly, score_from_labels, handing it both the truth and the labels — and check the result equals 2/3 (within floating-point tolerance, hence the .abs() < 1e-9).

The sibling test label_convenient_form_runs_inference_and_delegates then verifies the other path: the convenient form runs the model once to get labels, then delegates to the same pure form — and the two scores match. Together, the two tests pin down that the pure form never re-infers and that both paths agree. All of this lives in a companion folder beside the scoring implementation, exactly as the hygiene rules demand.

What both walk-throughs show

Neither of these tests is in some far-away tests/ directory. Each sits beside the code it tests, in a companion module declared at the end of the relevant mod.rs. Because of the tests-at-end rule, the implementation files read top-down, and because of the pure-dispatcher rule, their mod.rs files are tiny catalogues. When a future engineer wonders “is prefetching really overlapped?”, the test is one folder over.

6. Look inside

Here is the folder-module tree for a couple of the real modules we studied, showing how mod.rs sits as the dispatcher above implementation files and companion tests:

flowchart TD
    accTitle: Two folder modules showing the pure-dispatcher pattern
    accDescr: Each folder module has a mod.rs that acts as a pure dispatcher, implementation files for each concern, and companion test files. The kernels folder holds elementwise, grain, reductions and scaling; the execution folder holds context, modes, plan, resolve and the streaming driver.
    subgraph K["sciencekit_math :: kernels/ (folder module)"]
        direction TB
        K_RS["mod.rs<br/><i>pure dispatcher</i>"]
        K_ELEM["elementwise.rs"]
        K_GRAIN["grain.rs"]
        K_RED["reductions.rs"]
        K_SCAL["scaling.rs"]
        K_TESTS["kernels_tests.rs + grain_tests.rs<br/><i>companion tests</i>"]
    end

    subgraph E["sciencekit_common :: execution/ (folder module)"]
        direction TB
        E_RS["mod.rs<br/><i>pure dispatcher</i>"]
        E_CTX["context.rs"]
        E_MODES["modes.rs"]
        E_PLAN["plan.rs"]
        E_RES["resolve.rs"]
        E_DRV["streaming_driver.rs"]
        E_TESTS["execution_tests.rs + streaming_driver_tests.rs<br/><i>companion tests</i>"]
    end

    K_RS --> K_ELEM
    K_RS --> K_GRAIN
    K_RS --> K_RED
    K_RS --> K_SCAL
    K_RS -.-> K_TESTS

    E_RS --> E_CTX
    E_RS --> E_MODES
    E_RS --> E_PLAN
    E_RS --> E_RES
    E_RS --> E_DRV
    E_RS -.-> E_TESTS

    style K_RS fill:#ede9fe,stroke:#7c3aed,color:#3b0764,stroke-width:2px
    style E_RS fill:#ede9fe,stroke:#7c3aed,color:#3b0764,stroke-width:2px
    style K_TESTS fill:#ffe4e6,stroke:#e11d48,color:#7f1d1d,stroke-dasharray:5 5
    style E_TESTS fill:#ffe4e6,stroke:#e11d48,color:#7f1d1d,stroke-dasharray:5 5

The solid arrows show mod.rs declaring its implementation submodules; the dashed arrows show it declaring the #[cfg(test)] companion test modules. Notice that every solid child is a real implementation file and every dashed child is a test file — and mod.rs itself, at the top of each folder, holds no logic at all.

7. Recap

• A folder module (mod.rs inside a directory) is how Rust lets one conceptual module span multiple files once a single .rs file passes ~200 lines.
• mod.rs must be a pure dispatcher: it only mod-declares submodules and pub use-re-exports them. It carries no implementation logic, so it reads as a one-screen catalogue.
• Tests live in companion *_tests.rs modules beside the implementation — never in a distant global tests/ directory — and are declared with #[cfg(test)] at the end of mod.rs.
• Tests-at-end: any inline #[cfg(test)] block sits at the very end of its file, so the implementation reads top-down.
• This was a pure refactor: no behavior changed, every test stayed green. The value is purely in how the codebase is organized and how easy it is to navigate.

That is the discipline of module hygiene. In the next chapters we put this tidy filing cabinet to work as we climb back up to the algorithms — now that every folder has a clean label, we can always find where the math actually lives.

Parallel kernels and the streaming executor

Modern machines have many cores, but your data comes in all sizes. A library that blindly parallelizes everything wastes time on a 4-element array, and a library that never parallelizes ignores most of the hardware. And when a dataset is too large to hold in memory, the same library must keep the disk reader and the CPU working at the same time, not one after the other. This chapter explains how sciencekit decides how much parallelism each operation deserves, how the math kernels switch between sequential and parallel execution at the push of a single number, and how a streaming executor overlaps reading with computing so big data never idles the CPU.

1. The problem we were solving

Imagine a kitchen with several cooks. Each cook can chop vegetables independently — that is parallelism. But handing a single carrot to all four cooks is absurd: the time they spend arguing over who chops which slice is longer than just chopping it. The cooks are only worth using when there is a pile of carrots. Now add an oven: while one cook is waiting for the oven, the other cook can keep chopping the next tray. That overlap — preparing the next batch while the current one bakes — is what a streaming library must do when the vegetable pile is bigger than the counter can hold.

Before this change, sciencekit’s kitchen was incoherent. The library already had a resolver that produced an SKExecutionPlan with a parallelism field — how many threads to use — but nothing in the math layer actually read that field. The scaling kernel parallelized even a 4-element array (four cooks, one carrot), while the elementwise, binary, and reduction kernels never parallelized at all, no matter how big the data was. Parallelism was decided by accident, not by a plan.

Three things were missing. First, a decision rule: given the data size, the number of cores, and the known cost of each kernel, how many threads should an operation use — or should it stay sequential? Second, kernels that honor the plan: each math primitive must actually branch on the resolved parallelism. Third, a streaming executor: a reusable piece that runs the reader and the compute concurrently through a double buffer, so out-of-core data is processed without the CPU ever sitting idle waiting for the disk.

2. The decisions — and the roads not taken

3. The concepts, taught from zero

Rust concept 1: rayon and parallel iterators

rayon is a Rust library for easy data parallelism. Instead of writing a loop that visits each element, you hand rayon an iterator and it splits the work across a pool of threads, giving you the result as if you’d run the loop sequentially — but faster.

The clever part is that a normal iterator and a parallel iterator have the same shape, so switching costs one word. An iterator produces items one at a time via a method that yields an Option (an optional value: Some(item) or None). A parallel iterator (into_par_iter, par_iter) distributes those items across threads. The whole point of this chapter’s design is that the kernels use the same higher-order iteration for both paths, only changing the kind of iterator based on the resolved parallelism.

Rust concept 2: azip! and par_azip! — looping without index arithmetic

ndarray is the array library sciencekit uses. The macros azip! and par_azip! are its idiomatic way to visit all elements of one or more arrays together — the “apply f pairwise across these arrays” operation. azip! is sequential; par_azip! is the rayon-parallel version. They share the same syntax, which is exactly why the kernels can swap them with a one-line branch.

#![allow(unused)]
fn main() {
use ndarray::{Array, ArrayView, azip, par_azip};

let mut output = Array::zeros(input.dim());
if parallelism > 1 {
    par_azip!((out in &mut output, a in input) { *out = transform(*a); });
} else {
    azip!((out in &mut output, a in input) { *out = transform(*a); });
}
}

Inside the macro body, out is a mutable reference to each output element and a is a reference to each input element. The macro visits every element position once, pairing the output slot with its input element. This is a higher-order form — you provide the closure (the { ... } body) and the library provides the traversal. The project’s rule is never a manual index loop: writing output[i] = f(input[i]) by hand is forbidden, because the macro form is clearer, safer, and lets the same body run parallel with a one-word change.

Rust concept 3: closures capturing by move

A closure is an anonymous function you can define inline and pass around. |x| 2.0 * x + 1.0 takes one argument and returns twice it plus one. A closure can capture variables from the surrounding scope — “remember” them. In sciencekit’s kernels the closure must be Sync + Send:

#![allow(unused)]
fn main() {
pub fn sk_elementwise_transform<F: SKFloat>(
    input: &ArrayView<F, ndarray::Ix1>,
    transform: impl Fn(F) -> F + Sync + Send,
    parallelism: usize,
) -> Array<F, ndarray::Ix1>
}

• Sync means the value is safe to be referenced from multiple threads at once.
• Send means it is safe to be moved to another thread.

These bounds are required because par_azip! takes the closure and runs it on the rayon pool — across many threads simultaneously. The compiler enforces this: if you try to capture a value that is not Sync + Send, the code will not compile. This is Rust pushing the mistake out of the running program and into the compiler: a closure that cannot safely run in parallel is rejected before it ever runs.

Numerical concept 1: what “grain” means

Parallelism has a fixed dispatch overhead: splitting work, coordinating threads, and rejoining them all cost time. A single thread on a tiny array is therefore faster than many threads, because the overhead swamps the tiny amount of work. The grain is the minimum number of elements per thread that makes a thread worth spawning.

Think of it as the “amortization” point. If a thread must process at least g elements to pay for the cost of spawning it, then for a unit of u elements you can afford at most ceil(u / g) threads. The resolver computes:

threads = min(cores, ceil(unit_elements / grain))

• If unit_elements < grain, then ceil(unit / grain) = 1, so it runs on one thread — sequential.
• If unit_elements is huge, ceil climbs toward the core count, where it is capped.

The grain is measured per kernel by a calibration benchmark, not guessed. grain.rs records the measured crossover for each kernel on a specific machine and date — that is called provenance (where the number came from) — so a magic constant is never unexplained. The same grain applies to in-memory and streaming work, because dispatch overhead and per-element cost do not change just because the data came from disk.

Numerical concept 2: the “work unit” — whole dataset vs one batch

The grain needs a unit to be applied to: the number of elements we are parallelizing over. This is where mode matters:

• In-memory: the work unit is the whole dataset. unit = dataset_elements.
• Streaming: compute happens one batch at a time, so the work unit is one batch. unit = batch_size_hint.
• Memory-mapped: an arbitrary shard can be touched, so unit is treated as “always large” and resolves to cores.

This is why the same parallelism field means “threads for the current operation’s unit”, decided once at resolution and carried in the plan down to each kernel.

Numerical concept 3: double buffering — the producer/consumer overlap

Reading from disk (or any slow source) and computing are two different speeds. If they run strictly one after the other, one of them is always idle: the disk waits while the CPU computes, then the CPU waits while the disk reads.

Double buffering fixes this with a producer/consumer relationship:

• A producer (the reader) fetches batches from the source.
• A consumer (the compute side) processes them.
• Between them sits a bounded buffer — a fixed-size queue of batches.

The producer and consumer run on different threads. While the consumer computes batch k, the producer reads batch k+1 into the buffer. When the consumer finishes batch k, batch k+1 is already waiting — no idle time. The buffer is bounded (it holds at most buffer_depth prefetched batches, plus the one being computed) so the reader never runs unboundedly ahead.

The name “double buffering” refers to the simplest bound: buffer_depth = 1, so at most one batch is computing and one is prefetched — two resident batches. That is why the resolver checks $2 \times \text{batch} \le \text{available memory}$ before it will even allow a streaming plan: two batches must fit in RAM. The plan field is literally called buffer_depth, and 1 is double buffering.

Numerical concept 4: closures for per-element math, not shared accumulators

When a kernel reduces many rows into one answer (a sum), the naive parallel temptation is a single shared accumulator every thread adds to. That needs locking (slow) or atomics (also slow, and order-dependent). The design instead gives each chunk its own partial accumulator, computed in parallel with no sharing, then combines the partials at the end. This is lock-free and race-free by construction: no two threads ever touch the same memory while running. (The trade-off is that floating-point addition is not associative, so a parallel sum may differ from a sequential one in the last bits — which is why the tests compare with a tolerance, never bit-for-bit equality.)

4. Each object, explained

parallelism_for (private, but the heart of the decision)

What it is: A private helper inside the resolver (resolve.rs) that computes the parallelism for a resolved mode: min(cores, ceil(unit / grain)). It returns 1 if there are at most one core, if the unit is unknown or zero, or if the unit is below one grain.

Why it exists: It is the single, pure place where “how many threads?” is answered from the size of the work unit — whole dataset for in-memory, one batch for streaming, ∞ (so “always cores”) for memory-mapped shards. Because it is a pure function of the context, it is testable by injecting simulated values with no real machine involved (see grain_governs_resolution_across_regimes below).

What we rejected: A fixed global threshold, reading sysinfo inside the resolver, and deciding parallelism per-kernel at the point of execution. All of those either couple the decision to the physical machine or duplicate the policy. The number, the unit, and the grain all flow through this one formula, deterministically.

sk_elementwise_transform

What it is: The free-scope kernel sk_elementwise_transform(input, transform, parallelism) -> Array, which maps x → transform(x) over a 1-D vector and returns a brand-new owned array. It branches on parallelism: 1 runs the sequential azip!, > 1 dispatches par_azip!.

Why it exists: It is the most basic building block of “do something to every element.” It must exist in both sequential and parallel form because a transform closure can be cheap (e.g. multiply) or expensive (e.g. sin), so the right amount of parallelism changes. The iterator form (never a manual index loop) honors the project’s higher-order-function rule.

What we rejected: A single always-parallel version (a regression on the cheap-closure, small-array case) and two separate functions for the two paths (a duplicated API that pushes the decision onto every caller). One function, one parallelism number, two internal paths.

sk_binary_combine

What it is: The kernel sk_binary_combine(left, right, combine, parallelism) -> Array, pairing each element of left with the corresponding element of right and applying combine(left_value, right_value). It asserts the two inputs have equal length, seeds the output from left, then branches: > 1 uses a three-array par_azip!, 1 uses the sequential zip_mut_with.

Why it exists: Many operations (adding two matrices, elementwise products, comparisons) take two inputs and produce one output. It needs the same parallel/sequential switch as the unary transform. Its measured grain is the largest of the set (10⁷) precisely because a three-array par_azip! has the highest dispatch overhead — a measured fact recorded in grain.rs, not a guess.

What we rejected: Running combine on an index loop; writing into one input in place without seeding from left (that would destroy left); and a shared-accumulator approach. The three-array zip keeps all three buffers aligned by construction.

sk_axis_sum

What it is: The reduction kernel sk_axis_sum(input, axis, parallelism) -> Array, summing a 2-D array along one axis into a 1-D result. axis = 0 collapses rows (per-column sums, length = columns); axis = 1 collapses columns (per-row sums, length = rows). It has four paths: sequential and parallel for each axis.

Why it exists: Reductions are where parallelism is most valuable and most dangerous. The axis-0 parallel path uses per-chunk partial accumulation: rows are split into chunks, each chunk sums into its own partial vector, then the partials are combined. This never writes a shared accumulator concurrently. Axis-1 parallelizes trivially, one independent row per thread. The sequential fallbacks exist so small data (or parallelism = 1) avoids all dispatch.

What we rejected: Atomics on a shared accumulator (contention and nondeterministic order), and a single shared output array written by all threads (a data race). The partial-accumulation design is race-free by construction and only reorders the final floating-point additions by a few bits — which the tolerance-based assert_close contract in the tests accepts.

sk_run_streaming_driver

What it is: The free-scope executor sk_run_streaming_driver(source, state, update, plan) -> Result<(), SKError>. It spawns a dedicated producer thread that reads batches from an SKLazySource into a bounded double buffer, while the caller’s thread consumes each owned SKDataBatch and applies the algorithm-supplied update(batch, parallelism, state) closure. update returns an SKStreamDecision — Continue or Stop — to steer the pipeline.

Why it exists: It is the reusable embodiment of the I/O ∥ CPU overlap. Ordering, error propagation, final-batch delivery, and early stop are written and tested once; an algorithm gains streaming by writing only its update step. The compute side never blocks on input waits — the producer thread does the waiting, on its own thread, while rayon computes on the pool.

What we rejected: Reading and computing on the same thread (no overlap); letting the producer read unboundedly ahead (hence the bounded buffer_depth); and a Tokio-based pipeline for synchronous disk sources (an unnecessary dependency). A plain scoped thread is enough for the synchronous SKLazySource contract.

5. A real walk-through, using the tests

The most convincing way to learn is to read actual passing tests. Here are three: two from grain_tests.rs that exercise the decision rule, and one from kernels_tests.rs that proves a parallel reduction equals the sequential reference.

Test 1: grain_governs_resolution_across_regimes

This test proves that a unit below the grain stays sequential (parallelism = 1) and a unit far above it parallelizes — using only a simulated context, never the physical machine.

#![allow(unused)]
fn main() {
#[test]
fn grain_governs_resolution_across_regimes() {
    let grain = SK_AXIS_SUM_GRAIN;
    let below = context(Some((grain - 1) as u64), grain, 8);
    let above = context(Some((grain * 8) as u64), grain, 8);

    let below_plan = sk_resolve_execution_plan(SKExecutionMode::Automatic, &below).unwrap();
    let above_plan = sk_resolve_execution_plan(SKExecutionMode::Automatic, &above).unwrap();

    assert_eq!(below_plan.parallelism, 1);
    assert_eq!(above_plan.parallelism, 8);
}
}

Line by line:

1. let grain = SK_AXIS_SUM_GRAIN; — we grab the axis-sum kernel’s measured grain (100,000 elements) directly from the constant in grain.rs.
2. let below = context(Some((grain - 1) as u64), grain, 8); — the test helper context builds an SKExecutionContext by hand. The Some(grain - 1) is the dataset element count: one fewer than one grain. The grain argument is the kernel’s grain, and 8 is the simulated core count. So this context says: 99,999 elements, grain 100,000, 8 cores.
3. let above = context(Some((grain * 8) as u64), grain, 8); — a dataset of eight grains (800,000 elements) on the same 8 cores.
4. We resolve an Automatic plan for each. .unwrap() asserts resolution succeeds.
5. assert_eq!(below_plan.parallelism, 1) — because ceil(99,999 / 100,000) = 1, the plan is sequential even though 8 cores are available. The work is too small to pay for threads.
6. assert_eq!(above_plan.parallelism, 8) — ceil(800,000 / 100,000) = 8, clamped to the 8 cores, so it parallelizes fully.

The lesson: the grain is the throttle. The same resolver, same cores — only the size changed, and the parallelism followed the min(cores, ceil(unit / grain)) rule exactly. Nothing about a real machine was read.

Test 2: every_kernel_exposes_a_positive_grain

This short test guards the provenance: every kernel must have a real, positive grain constant — never a zero or a placeholder that would silently disable parallelization.

#![allow(unused)]
fn main() {
#[test]
fn every_kernel_exposes_a_positive_grain() {
    for grain in [
        SK_ELEMENTWISE_TRANSFORM_GRAIN,
        SK_BINARY_COMBINE_GRAIN,
        SK_AXIS_SUM_GRAIN,
        SK_SCALE_IN_PLACE_GRAIN,
    ] {
        assert!(grain > 0);
    }
}
}

Line by line:

1. The test iterates over the four measured grain constants from grain.rs.
2. assert!(grain > 0) — each must be strictly positive. A grain = 0 would make ceil(unit / 0) a division-by-zero hazard; a grain of usize::MAX would effectively disable parallelism forever. This test pins the constants to a sane, measured range.

The lesson: these are measured numbers with provenance (machine, date, protocol) recorded in grain.rs, and this test keeps them honest. A kernel whose grain was never calibrated cannot silently ship a wrong constant.

Test 3: parallel_axis_zero_sum_matches_sequential_reference

This test proves the race-free parallel axis-0 reduction equals the sequential column sums within floating-point tolerance — across many rows, where a naive shared accumulator could race.

#![allow(unused)]
fn main() {
#[test]
fn parallel_axis_zero_sum_matches_sequential_reference() {
    let (rows, cols) = (4096, 256);
    let input: Array2<f64> =
        Array2::from_shape_fn((rows, cols), |(i, j)| (i as f64) * 0.001 + (j as f64) * 0.5);
    let sequential = sk_axis_sum(&input.view(), 0, 1);
    let parallel = sk_axis_sum(&input.view(), 0, 8);
    assert_close(&sequential, &parallel);
}
}

Line by line:

1. let (rows, cols) = (4096, 256); — a wide, tall matrix: 4,096 rows and 256 columns.
2. Array2::from_shape_fn((rows, cols), |(i, j)| ...) — build the array by computing each cell from its position (i, j). This avoids listing 1M numbers by hand and gives every element a deterministic value. The closure |(i, j)| (i as f64) * 0.001 + (j as f64) * 0.5 blends row and column contributions so the columns are genuinely different.
3. let sequential = sk_axis_sum(&input.view(), 0, 1); — reduce along axis 0 (per-column sums) with parallelism = 1, which runs the plain sequential path.
4. let parallel = sk_axis_sum(&input.view(), 0, 8); — the same reduction with parallelism = 8, which takes the parallel path: rows split into chunks, each chunk sums into its own partial vector, then the partials are combined.
5. assert_close(&sequential, &parallel) — compare the two within a relative floating-point tolerance (1e-9, scaled by magnitude). Because the parallel path reorders the floating-point additions (chunk partials summed in a different order), the result may differ in the last bits — but never beyond tolerance. This is the contract the design documents: parallel is correct, not bit-identical.

The lesson: the parallel path produces the right column sums across 4,096 rows without any shared-accumulator race — and the tolerance-based comparison is what makes “parallel equals sequential” a testable guarantee rather than a hope.

6. Look inside

Here is the full pipeline, from resolution to a parallel kernel, and the double-buffer producer/consumer that keeps streaming from ever idling the CPU:

flowchart TD
    accTitle: The resolve, parallel kernel and streaming pipeline
    accDescr: Resolution turns an execution context into an execution plan with parallelism and batch size. Kernels branch on parallelism to use a parallel or sequential iterator. A streaming plan hands the parallelism to a double-buffered producer and consumer, so reading and computing overlap.
    subgraph resolve[Resolution — pure & deterministic]
        CTX[SKExecutionContext<br/>memory · cores · unit size · grain] --> RES[sk_resolve_execution_plan]
        RES --> PLAN[SKExecutionPlan<br/>mode · parallelism · batch_size · buffer_depth]
    end

    PLAN -->|parallelism & batch_size| KERNEL
    subgraph compute[Parallel kernels]
        KERNEL{parallelism > 1?}
        KERNEL -- yes --> PAR[par_azip! across the compute pool]
        KERNEL -- no --> SEQ[azip! / zip_mut_with — sequential]
        PAR --> OUT[owned result]
        SEQ --> OUT
    end

    PLAN -->|OutOfCoreStreaming| DRIVER
    subgraph stream[Streaming executor — producer/consumer]
        PROD[Producer thread · reads batch k+1] --> BUF[Double buffer<br/>buffer_depth = 1]
        BUF --> CON[Consumer · update batch k → state]
        CON -->|Continue / Stop| CON
    end

    style resolve fill:#fef3c7,stroke:#d97706,color:#78350f
    style compute fill:#f0fdf4,stroke:#16a34a,color:#14532d
    style stream fill:#e0f2fe,stroke:#0284c7,color:#0c4a6e

Reading the left branch: resolution decides parallelism from the context’s unit size and the kernel’s grain, and each kernel branches on that one number to pick its sequential or parallel iterator. Reading the right branch: a streaming plan hands the same parallelism to the driver, whose producer thread prefetches batch k+1 while the consumer computes batch k on the rayon pool — compute never blocks on input, and the reader never runs more than buffer_depth ahead.

Here is that double-buffer overlap as a timeline, showing the producer and consumer advancing in step:

sequenceDiagram
    accTitle: Double-buffered streaming timeline
    accDescr: The producer reads batches into a double buffer while the consumer computes, so reading batch one overlaps computing batch zero and the reader never runs more than buffer depth ahead.
    participant P as Producer (reader)
    participant B as Double buffer
    participant C as Consumer (compute)

    P->>B: enqueue batch 0
    B->>C: deliver batch 0
    Note over P,C: compute batch 0 overlaps read batch 1
    P->>B: enqueue batch 1
    C-->>C: update(batch 0)
    B->>C: deliver batch 1
    C-->>C: update(batch 1)
    Note over P,C: compute batch 1 overlaps read batch 2
    P->>B: enqueue batch 2
    B-->>C: ... continue

The producer and the consumer advance on different threads, so while the consumer runs update(batch k), the producer is already enqueuing batch k+1. That is the I/O ∥ CPU overlap the spec streaming-executor requires, bounded by the plan’s buffer_depth.

7. Recap

• Parallelism is a decision, not a kernel trait. sk_resolve_execution_plan computes min(cores, ceil(unit / grain)) once, and every kernel honors the resulting parallelism with a single azip!/par_azip! branch.
• The grain is measured, with provenance. Each kernel’s grain (minimum elements per thread) comes from a calibration benchmark on a specific machine and date — not a guessed constant — and a positive-grain test keeps the constants honest.
• The work unit changes with the mode. In-memory parallelizes over the whole dataset, streaming over one batch, memory-mapped over arbitrary shards — but the same grain rule governs all three.
• Reductions are parallelized race-free. Axis-0 sums use per-chunk partial accumulation, never a shared accumulator, and compare against the sequential reference within float tolerance rather than bit-identity.
• The streaming executor overlaps I/O and compute. sk_run_streaming_driver runs a producer thread into a bounded double buffer while the consumer computes on the rayon pool, delivering ordered, owned batches and stopping early on request.

Next chapter: now that an operation can parallelize its math and stream its data, we look at how the whole workspace is shaped into crates — the boundaries that keep kernels, execution, and the umbrella library decoupled and testable.

Maintaining this book

A textbook is only worth reading if it stays true. This book explains why sciencekit is built the way it is, and sciencekit is not finished — the roadmap (docs/PRD.md) maps out new algorithms in phases that stretch from preprocessing all the way to trees, ensembles, clustering and Python bindings. Every one of those algorithms will change the library, and every change risks turning a sentence in this book into a lie. This chapter is the contract that keeps the book truthful: a small set of rules that decide when to add a chapter, when to edit one, and how every update flows safely from a branch to the published site.

1. The problem we were solving

Documentation rots. It is not a bug in any single sentence; it is a slow drift. Someone adds a new algorithm, nobody updates the matching page, and six months later the book describes a SKKNeighborsClassifier that behaves nothing like the code. A reader who trusts the book builds the wrong mental model, and the wrong mental model produces subtle, expensive bugs.

You cannot stop change, so the only fix is to make keeping up a mechanical habit. sciencekit already commits to this habit at the repository level — its AGENTS.md requires that every code change ships a companion documentation change on its own branch, merged separately. This chapter is where that general rule becomes a concrete procedure for this book: when a feature lands you add a chapter; when behavior changes you update the chapter; and both are done through the same worktree-and-branch machinery that keeps main clean.

2. The decisions — and the roads not taken

3. The concepts, taught from zero

What “shipping a doc change” means

In most projects you might finish a feature and also remember to update the docs. In sciencekit that is not optional and not deferred: AGENTS.md says every code change is paired with a documentation change that updates the mdBook for the newly introduced behavior, “on its own doc branch and its own PR merged into docs/documentations”. The doc branch is created at the same time as the code change, not after it.

The two documentation branches

• docs/documentations is the landing branch: the finished book lives here, and the GitHub Actions workflow (.github/workflows/deploy-documentation.yml) builds and deploys it to GitHub Pages on every push.
• docs/<name> is a working branch. You never write on docs/documentations directly; you open a worktree for your change and merge it in when it is ready.

This is deliberately separate from the -openspec planning branch, which only holds the OpenSpec planning artifacts. Documentation (the book) and planning (the spec) are two different concerns and use two different branches.

What “worktree” means here

A git worktree lets you check out a branch into a separate folder without disturbing your current checkout. Every docs change lives in temporary/worktrees/docs/<name> on branch docs/<name>. temporary/ is git-ignored, so these working copies are never committed — they are scratch space that keeps each change isolated.

4. Each object, explained

The chapter file

What it is — a markdown file under docs/src/development-explained/, one per chapter, following the anatomy in the chapter template (CHAPTER_TEMPLATE.md).

Why it exists — each chapter is the durable record of one concept’s design rationale, taught from zero for a reader who knows programming but not Rust or numerical computing.

What we rejected — putting everything in one huge page; a single page cannot teach or be maintained.

The worktree and branch

What it is — temporary/worktrees/docs/<name> on branch docs/<name>.

Why it exists — it isolates a doc change so the deployable docs/documentations branch only ever receives finished, reviewed prose.

What we rejected — editing docs/documentations directly; unreviewed edits reach the published site immediately.

The build gate

What it is — mdbook build docs, run locally before pushing.

Why it exists — it validates every chapter and its links before the change is merged; the book’s config uses create-missing = false, so a broken link fails the build rather than silently generating an empty page. The gate must pass with zero warnings.

What we rejected — relying on the remote deploy workflow to catch errors; that would let broken content reach the published book.

5. A real walk-through: adding an algorithm’s chapter

Here is the exact procedure, using a hypothetical new algorithm from a later roadmap phase (say a tree model arriving in Phase 4).

Step 1 — open the worktree. Create the docs change alongside the code change:

git worktree add temporary/worktrees/docs/sk-decision-tree docs/sk-decision-tree

Step 2 — write the chapter. Create docs/src/development-explained/chNN-*.md following the template: the plain-language problem, the decisions table, concepts taught from zero, each object explained, and a walk-through built from the real *_tests.rs tests. Keep code identifiers exactly as they are in the source.

Step 3 — register the chapter. Add it to the table of contents in docs/src/SUMMARY.md under the right part of the book, in the position that matches the reading order.

Step 4 — run the gate.

mdbook build docs

Fix every warning. A warning here means the published book could be broken.

Step 5 — merge. Open a PR for branch docs/sk-decision-tree into docs/documentations. When it merges, the deploy workflow rebuilds and publishes the book automatically.

Updating, not forking

If the change is a bug fix or refactor, replace the affected part of the existing chapter — the decisions table, the object descriptions, the walk-through — and adjust any prose that now contradicts the code. Do not copy the chapter and add a “v2”. The book stays truthful because there is exactly one place where a given behavior is explained, and it is always the current explanation.

6. Look inside

flowchart LR
    accTitle: How a code change flows into the documentation
    accDescr: A code change opens a code worktree branch and a docs worktree branch. The code branch becomes a code pull request. The docs branch must pass mdbook build with zero warnings, then becomes a doc pull request merged into docs/documentations, which deploys to GitHub Pages.
    C[code change] --> B1[worktree code branch]
    C --> B2["worktree docs branch"]
    B1 --> P1[code PR]
    B2 --> G["mdbook build docs<br/>zero warnings"]
    G --> P2["doc PR → docs/documentations"]
    P2 --> D["deploy → GitHub Pages"]

The docs branch is created at the same moment as the code branch, runs through its own build gate, and lands on docs/documentations through its own PR — the same “two coexisting worktrees, two separate PRs” discipline that governs code.

7. Recap

• A new feature or algorithm → add a new chapter, registered in SUMMARY.md.
• A bug fix or refactor → update the existing chapter in place; never add a parallel twin.
• Every code change ships its doc change on a docs/<name> branch merged into docs/documentations — created alongside the code change, never deferred.
• mdbook build docs must pass with zero warnings before you push; broken links fail the build, and the deploy happens automatically on merge.
• Keep counts, phases and status in sync with docs/PRD.md.

Next chapter: we step down from the meta and back into the code, and see how one of those future algorithms will actually be built.