Metadata-Version: 2.4
Name: biomimetic-framework
Version: 0.2.0
Summary: Biomimetic Lookup Cell — parametric geodesic lookup framework for machine learning
License: MIT
Project-URL: Homepage, https://github.com/SariItani/bmc
Project-URL: Documentation, https://SariItani.github.io/bmc
Project-URL: Repository, https://github.com/SariItani/bmc
Project-URL: Issues, https://github.com/SariItani/bmc/issues
Keywords: machine-learning,lookup-table,voronoi,geodesic,gradient-free
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: Topic :: Scientific/Engineering :: Artificial Intelligence
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.9
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Requires-Python: >=3.9
Description-Content-Type: text/markdown
Requires-Dist: numpy>=1.24
Requires-Dist: scipy>=1.10
Provides-Extra: image
Requires-Dist: Pillow>=9.0; extra == "image"
Requires-Dist: torch>=2.0; extra == "image"
Requires-Dist: opencv-python>=4.7; extra == "image"
Provides-Extra: gpu
Requires-Dist: torch>=2.0; extra == "gpu"
Provides-Extra: viz
Requires-Dist: rich>=12.0; extra == "viz"
Requires-Dist: graphviz>=0.20; extra == "viz"
Provides-Extra: dev
Requires-Dist: pytest>=7.0; extra == "dev"
Requires-Dist: pytest-cov; extra == "dev"

# BMC Framework — API Reference

**Version 0.1.13**

---

## Overview

The `bmc` library is a gradient-free machine learning framework built around the
**Biomimetic Lookup Cell** primitive. Instead of learned weights, it stores empirical
probability distributions in Voronoi cells ordered along a geodesic path through the
data manifold.

### Top-level exports

All primary classes and functions are importable directly from `bmc`:

```python
# High-level API
from bmc import BMCClassifier, BMCRegressor, save, load

# Composition
from bmc import BMCLayer, BMCNetwork, SequentialBMC, BMCGate, BMCSkipConnection, SequentialRVQ

# Core algorithms
from bmc import fit_rvq, reconstruct_rvq

# Primitive
from bmc import BMCCell

# Observation spaces
from bmc import (
    ImageObservationSpace, TabularObservationSpace,
    TimeSeriesObservationSpace, TextObservationSpace,
)

# Metrics
from bmc import L2Metric, FrobeniusMetric, CosineMetric, SpectralMetric, NuclearMetric

# Interpolators
from bmc import LinearInterpolator, SplineInterpolator, KernelInterpolator, RDDAInterpolator

# FPGA export
from bmc import export_for_fpga

# Visualization
from bmc.viz import summary, diagram, scatter, architecture, trace, parametric, riemannian
```

### Class hierarchy

```
ObservationSpace          — decomposes raw input into zone feature vectors
    ImageObservationSpace
    TabularObservationSpace
    TimeSeriesObservationSpace
    TextObservationSpace

SimilarityMetric           — computes distances between feature vectors
    L2Metric
    FrobeniusMetric
    CosineMetric
    SpectralMetric
    NuclearMetric

Interpolator (ABC)         — geodesic interpolation strategy
    LinearInterpolator     — piecewise-linear (fast baseline)
    SplineInterpolator     — cubic spline (C2 continuous)
    KernelInterpolator     — IDW or RBF kernel
RDDAInterpolator           — feature-space Delta-Adjust (not an Interpolator subclass)

BMCCell                    — single-zone Voronoi lookup unit (the atomic primitive)

BMCLayer                   — N BMCCells in parallel over N zones of one input

BMCNetwork                 — layers stacked sequentially, arbitrary depth
    SequentialBMC          — convenience constructor for BMCLayer → BMCCell

BMCGate                    — geodesic-position router (regression BMCCell on arc-length)

BMCSkipConnection          — adaptive skip connection for representation collapse
                             L1 + BMCGate + specialist L2 (entropy-weighted seeding)

SequentialRVQ              — sequential BMC via Residual Vector Quantization
                             auto-determines depth, classification + regression

fit_rvq(cells, X)          — core RVQ training algorithm (trains cells on residuals)
reconstruct_rvq(cells, X)  — accumulated RVQ reconstruction

BMCClassifier              — sklearn-compatible high-level classifier
BMCRegressor               — sklearn-compatible high-level regressor

save(model, path)          — serialize a fitted model to .bmc file
load(path, ...)            — restore a fitted model from .bmc file

export_for_fpga(cell, dir) — export a fitted BMCCell to FPGA Verilog + ROM hex

summary(model)             — CLI architecture summary (rich tree)
diagram(model, path)       — Graphviz architecture diagram (SVG/PDF)
scatter(model, X, y, path) — UMAP/PCA scatter with Voronoi + geodesic
architecture(model, X, y)  — composite: architecture with embedded scatter plots
trace(model, X, y, query)  — prediction trace with highlighted cells + geodesic edge
parametric(cell, path)     — geodesic curvature plot with class color gradient
riemannian(cell, path)     — 3D Riemannian manifold in joint (D+C) space
```

### Data flow

```
raw input
    └─► ObservationSpace.extract()  →  (N_zones, M, D)  zone feature vectors
            └─► BMCLayer.fit()      →  N_zones × BMCCell, each with K seeds + LUT
                    └─► BMCLayer.transform()  →  (M, N_zones, C)  vote matrix
                            └─► BMCNetwork    →  (M, C)  refined distribution
                                    └─► argmax  →  (M,)  predictions
```

### Dtype policy (since 0.1.12)

Output dtype follows input dtype throughout the pipeline. **All float types
are supported**: float16, float32, float64 (and any other `np.floating`
subtype). Non-float inputs (int8, uint8, int32, etc.) are cast to float32.

- float16 input → float16 seeds, LUT, arc-lengths, vote matrices, predictions
- float32 input → float32 throughout
- float64 input → float64 throughout
- int/uint input → cast to float32, then float32 throughout

Internal computation may use higher precision where needed (e.g., Lloyd's
centroid accumulator), but results are cast back to the input dtype.

Where this document annotates return types as `float64` or `float32`,
the actual dtype follows the input unless otherwise noted.

---

## bmc.core.cell — BMCCell

```python
class bmc.core.cell.BMCCell(
    n_clusters,
    metric=None,
    interpolate=False,
    interpolator=None,
    rng=None,
    batch_k=None,
)
```

The atomic primitive of the BMC framework. Encapsulates a single zone's complete
pipeline: manifold seeding → geodesic ordering → LUT construction → inference.

Analogous to `nn.Linear` in PyTorch — the fundamental building block from which
all higher-level structures are composed.

**Parameters**

| Name | Type | Default | Description |
|---|---|---|---|
| `n_clusters` | `int` | required | Number of Voronoi cells K. Acts as the memory budget for this zone. Each cell stores one row of the LUT — an empirical probability distribution `P(class \| cell)`. |
| `metric` | `SimilarityMetric` | `L2Metric()` | Distance metric for seed selection and Voronoi assignment. |
| `interpolate` | `bool` | `False` | Backward-compatible flag. If `True` and no `interpolator` is provided, uses `SplineInterpolator`. Ignored if `interpolator` is explicitly set. |
| `interpolator` | `Interpolator` or `RDDAInterpolator` | `None` | Interpolation strategy. Pass `None` for hard lookup, or an interpolator instance. Overrides `interpolate` flag. |
| `rng` | `np.random.Generator` | `default_rng(42)` | Random generator for kmeans++ initialization. Pass an explicit generator for full reproducibility. |
| `batch_k` | `int` or `None` | `None` | Seed-batch size for pairwise distance computation during Voronoi assignment. `None` auto-computes from 75% of available memory (system RAM for CPU metrics, GPU VRAM for CUDA metrics). User override accepted; values exceeding 99% of available memory are clamped with a warning. Passed through to `kmeanspp_init` and `assign_clusters` during `fit()`. |

**Attributes** *(set after `fit()`)*

| Name | Type | Description |
|---|---|---|
| `seeds_` | `ndarray, (K, D)` | Fitted seed vectors in TSP geodesic order. |
| `lut_` | `ndarray, (K, C)` | Lookup table in TSP order. `lut_[k]` is a probability distribution over C classes summing to 1.0. |
| `tsp_path_` | `ndarray, (K,)` | Integer permutation giving the TSP ordering of the original seeds. |
| `arc_lengths_` | `ndarray, (K,)` | Normalized arc-length positions of seeds along the geodesic. |
| `n_classes_` | `int` | Number of output classes. |
| `is_fitted_` | `bool` | `True` after `fit()` has been called. |

---

### BMCCell.fit

```python
BMCCell.fit(vectors, labels, n_classes=None, sample_weights=None) → BMCCell
```

Fit the cell on training vectors and labels.

**Steps executed:**
1. Initialize K seeds via kmeans++ from `vectors` (optionally weighted by `sample_weights`)
2. Refine seeds by Lloyd's iteration (assign → recompute means → repeat until convergence)
3. Assign each training vector to its nearest converged seed (Voronoi assignment)
4. Build LUT: for each cluster k, compute `P(class | k)` from training counts
5. Compute TSP geodesic ordering of seeds
6. Reorder `seeds_` and `lut_` by TSP path
7. Compute arc-length positions and fit interpolator (if set)

**Parameters**

| Name | Type | Description |
|---|---|---|
| `vectors` | `ndarray, (M, D)` | Training feature vectors for this zone. float32 or float64; non-float types are cast to float32. Output dtype (seeds, LUT) follows input dtype. |
| `labels` | `ndarray, (M,), int` | Class labels. Values must be in `[0, n_classes)`. |
| `n_classes` | `int`, optional | Total number of classes. If `None`, inferred as `labels.max() + 1`. |
| `sample_weights` | `ndarray, (M,)`, optional | Per-point importance weights for seed initialization. Biases seed placement toward high-weight regions (e.g., high-entropy cells from a previous layer). `None` → uniform weights (standard kmeans++). |

**Returns** `self` — enables chaining: `cell.fit(X, y).predict(X_test)`

**Notes**

- Empty Voronoi cells (no training samples assigned) receive a uniform distribution.
- Seeding uses kmeans++ initialization followed by Lloyd's iterative refinement until convergence. No gradients are used — the optimization is purely geometric (mean recomputation).
- When M > 10×K, Lloyd's iterations use **mini-batch sampling** (`lloyd_batch = min(M, 10*K)`) for efficiency. A final full-dataset pass ensures exact centroids. Accuracy impact is negligible (≤0.04% on overlapping 10-class data).

---

### BMCCell.transform

```python
BMCCell.transform(vectors) → ndarray, (M, C)
```

Map input vectors to probability distributions over classes.

If `interpolate=False` (default): returns `lut_[nearest_seed(vector)]` — a hard
lookup of the nearest seed's stored distribution.

If `interpolate=True`: selects the closer of the two TSP neighbors (symmetric
selection), projects the query onto the geodesic segment between the nearest seed
and the selected neighbor, then blends the two distributions proportionally to the
projection parameter t.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `vectors` | `ndarray, (M, D)` | Input feature vectors. |

**Returns**

| Name | Type | Description |
|---|---|---|
| `distributions` | `ndarray, (M, C)` | Per-sample probability distributions. Each row sums to 1.0. Dtype follows input. |

**Raises** `RuntimeError` if called before `fit()`.

---

### BMCCell.predict

```python
BMCCell.predict(vectors) → ndarray, (M,)
```

Predict class labels as `argmax` of `transform(vectors)`.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `vectors` | `ndarray, (M, D)` | Input feature vectors. |

**Returns**

| Name | Type | Description |
|---|---|---|
| `predictions` | `ndarray, (M,), int32` | Predicted class index per sample. |

**Raises** `RuntimeError` if called before `fit()`.

---

### BMCCell.fit_regression

```python
BMCCell.fit_regression(vectors, targets, sample_weights=None) → BMCCell
```

Fit the cell for regression on continuous targets. Identical pipeline to `fit()`
except the LUT stores mean target values per cell `(K, 1)` instead of class
probability distributions `(K, C)`. Rows are NOT normalized.

Accepts `sample_weights` for weighted seed initialization, same as `fit()`.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `vectors` | `ndarray, (M, D)` | Training feature vectors. Dtype policy applies. |
| `targets` | `ndarray, (M,), float64` | Continuous target values. |
| `sample_weights` | `ndarray, (M,)`, optional | Per-point importance weights for seed initialization. Same as `fit()`. |

**Returns** `self`

---

### BMCCell.predict_regression

```python
BMCCell.predict_regression(vectors) → ndarray, (M,)
```

Predict continuous values (regression mode). Uses hard lookup or geodesic
interpolation if an interpolator is fitted.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `vectors` | `ndarray, (M, D)` | Input feature vectors. |

**Returns** `ndarray, (M,)` — continuous predicted values. Dtype follows input.

**Raises** `RuntimeError` if called before `fit_regression()`.

---

### BMCCell.trace_query

```python
BMCCell.trace_query(query) → dict
```

Trace a single query through the cell's geodesic. Returns the nearest seed
index, chosen geodesic neighbor, parametric position t on the edge, and
arc-length position.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `query` | `ndarray, (D,)` | Single input vector. |

**Returns** `dict` with keys: `i` (nearest seed), `j` (geodesic neighbor),
`t` (parametric position 0-1), `arc` (global arc-length), `dist_i`, `dist_j`.

**Raises** `RuntimeError` if called before `fit()`.

---

### Aliases

`BMCCell.fit_classification` is an alias for `BMCCell.fit`.
`BMCCell.predict_classification` is an alias for `BMCCell.predict`.

---

### BMCCell.purity

```python
BMCCell.purity() → ndarray, (K,)
```

Return the purity of each Voronoi cell.

Purity = `max(lut_[k])` — the probability mass on the dominant class.
A pure cell has purity 1.0 (all assigned training samples belong to one class).
A maximally mixed cell has purity `1/C`.

**Returns**

| Name | Type | Description |
|---|---|---|
| `purities` | `ndarray, (K,)` | Per-cell purity values in `[1/C, 1.0]`. Dtype follows input. |

**Raises** `RuntimeError` if called before `fit()`.

---

### BMCCell.dominant_classes

```python
BMCCell.dominant_classes() → ndarray, (K,)
```

Return the dominant (most probable) class for each Voronoi cell.

**Returns**

| Name | Type | Description |
|---|---|---|
| `dominant` | `ndarray, (K,), int32` | `dominant[k] = argmax(lut_[k])`. |

**Raises** `RuntimeError` if called before `fit()`.

---

### BMCCell.get_state / set_state

```python
BMCCell.get_state() → dict
BMCCell.set_state(state) → BMCCell
```

Serialize and restore cell state. Used internally by `save()` / `load()`.
The state dict contains all fitted arrays and hyperparameters. No pickle — arrays
are stored as numpy arrays, scalars as JSON-native types.

---

## bmc.composition.layer — BMCLayer

```python
class bmc.composition.layer.BMCLayer(
    observation_space,
    n_clusters,
    metric=None,
    interpolate=False,
    rng=None,
    batch_k=None,
)
```

Parallel composition of N_zones `BMCCell` instances, one per zone of the input.

Each cell independently observes one zone, builds its own LUT, and produces a
probability distribution. The distributions are accumulated by soft voting (summed
then argmax-ed) to produce the final prediction.

**Parameters**

| Name | Type | Default | Description |
|---|---|---|---|
| `observation_space` | `ObservationSpace` | required | Defines how raw inputs are decomposed into zone feature vectors. |
| `n_clusters` | `int` | required | Number of Voronoi cells per zone (K). All zones share the same K. |
| `metric` | `SimilarityMetric` | `L2Metric()` | Distance metric applied uniformly across all zones. |
| `interpolate` | `bool` | `False` | Passed to each `BMCCell`. If `True`, all zones use geodesic interpolation. |
| `rng` | `np.random.Generator` | `default_rng(42)` | Each zone's cell receives a distinct derived seed from this generator. |
| `batch_k` | `int` or `None` | `None` | Seed-batch size. When `None` and parallel fitting is used, `BMCLayer.fit()` precomputes a safe value by dividing 75% of available memory (RAM or VRAM) by the number of worker threads, preventing concurrent overcommit. The computed value is propagated to all zone `BMCCell` instances. |

**Attributes** *(set after `fit()`)*

| Name | Type | Description |
|---|---|---|
| `layers_` | `list of BMCCell` | Fitted cells, one per zone. `layers_[z]` is the cell for zone z. |
| `n_classes_` | `int` | Number of output classes. |
| `is_fitted_` | `bool` | `True` after `fit()`. |

---

### BMCLayer.fit

```python
BMCLayer.fit(inputs, labels, n_classes=None) → BMCLayer
```

Fit all zone cells on training data.

Zone fitting is parallelized across CPU cores via `ThreadPoolExecutor` when
there are 4 or more zones. Each zone's `BMCCell` is fitted independently on
its own zone vectors. Threads share memory (no data duplication), and NumPy
releases the GIL during compute-heavy operations for real parallelism.
Falls back to sequential fitting for fewer than 4 zones.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `inputs` | `ndarray` | Raw training inputs. Shape depends on the `ObservationSpace`. For images: `(M, H, W)`. For tabular: `(M, F)`. |
| `labels` | `ndarray, (M,), int` | Class labels. |
| `n_classes` | `int`, optional | Total number of classes. Inferred from `labels` if `None`. |

**Returns** `self`

---

### BMCLayer.transform

```python
BMCLayer.transform(inputs) → ndarray, (M, N_zones, C)
```

Return per-zone probability distributions — the **vote matrix**.

The vote matrix is the intermediate representation that feeds into a second-layer
`BMCCell` or `BMCNetwork`. Each zone independently votes on the class; the vote
matrix records those votes before aggregation.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `inputs` | `ndarray` | Raw inputs in the same format as `fit()`. |

**Returns**

| Name | Type | Description |
|---|---|---|
| `vote_matrix` | `ndarray, (M, N_zones, C)` | `vote_matrix[i, z, c]` = probability assigned to class c by zone z for sample i. Each `vote_matrix[i, z, :]` sums to 1.0. float32. |

**Raises** `RuntimeError` if called before `fit()`.

---

### BMCLayer.predict

```python
BMCLayer.predict(inputs) → ndarray, (M,)
```

Predict class labels by summing votes across all zones and taking argmax.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `inputs` | `ndarray` | Raw inputs. |

**Returns**

| Name | Type | Description |
|---|---|---|
| `predictions` | `ndarray, (M,), int32` | Predicted class per sample. |

---

### BMCLayer.predict_proba

```python
BMCLayer.predict_proba(inputs) → ndarray, (M, C)
```

Return normalized probability distributions (accumulated votes, row-normalized).

**Returns**

| Name | Type | Description |
|---|---|---|
| `probabilities` | `ndarray, (M, C)` | Each row sums to 1.0. Dtype follows input. |

---

## bmc.composition.sequential — BMCNetwork

```python
class bmc.composition.sequential.BMCNetwork(layers)
```

General N-layer composable BMC network. Any layer implementing
`fit(X, y)` / `transform(X)` / `predict(X)` can be placed at any position.

Each layer's `transform()` output is flattened to `(M, D)` and passed as input
to the next layer. Labels are passed unchanged to every layer's `fit()` call —
every layer is supervised by the same ground-truth labels.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `layers` | `list` | Ordered list of layers. Minimum length 1. Each must implement `fit`, `transform`, and `predict`. |

**Raises** `ValueError` if `layers` is empty.
**Raises** `TypeError` if any layer is missing a required method.

**Attributes** *(set after `fit()`)*

| Name | Type | Description |
|---|---|---|
| `layers_` | `list` | Fitted layers in order. |
| `n_classes_` | `int` | Number of output classes. |
| `is_fitted_` | `bool` | `True` after `fit()`. |

**Notes on topology**

- `BMCLayer` at position 0 receives raw inputs; its `ObservationSpace` must match the input domain.
- `BMCLayer` at position > 0 receives a flat `(M, D)` array. Its `ObservationSpace` is automatically replaced with a `TabularObservationSpace` matching the input dimension.
- `BMCCell` at any position receives a flat `(M, D)` array.
- A layer that is already fitted when `fit()` is called will **not** be re-fitted.

**Common topologies**

```python
# 2-layer: parallel zones → single refinement cell
BMCNetwork([
    BMCLayer(obs, n_clusters=1024),
    BMCCell(n_clusters=256),
])

# 3-layer: parallel → refine → refine
BMCNetwork([
    BMCLayer(obs, n_clusters=1024),
    BMCCell(n_clusters=256),
    BMCCell(n_clusters=64),
])

# 3-layer: parallel → parallel (on vote matrix) → single
BMCNetwork([
    BMCLayer(obs_image, n_clusters=1024),
    BMCLayer(obs_any_n_zones_4, n_clusters=256),   # auto-rewrapped
    BMCCell(n_clusters=64),
])
```

---

### BMCNetwork.fit

```python
BMCNetwork.fit(X, y, n_classes=None) → BMCNetwork
```

Fit all layers sequentially.

Layer 0 is fitted on `X`. Layer 1 is fitted on `layer0.transform(X)` (flattened).
Layer k is fitted on `layer_{k-1}.transform(...)` (flattened). Labels `y` are
passed to every layer.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `X` | `ndarray` | Raw training input for layer 0. |
| `y` | `ndarray, (M,), int` | Class labels. |
| `n_classes` | `int`, optional | Inferred from `y` if `None`. |

**Returns** `self`

---

### BMCNetwork.transform

```python
BMCNetwork.transform(X) → ndarray, (M, C)
```

Run the full forward pass and return the final layer's output distribution.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `X` | `ndarray` | Raw input (same format as `fit()`). |

**Returns**

| Name | Type | Description |
|---|---|---|
| `output` | `ndarray, (M, C)` | Final layer's probability distributions. |

**Raises** `RuntimeError` if called before `fit()`.

---

### BMCNetwork.transform_at

```python
BMCNetwork.transform_at(X, depth) → ndarray
```

Run the forward pass up to and including layer `depth`.

Useful for inspecting intermediate representations — the gradient-free analogue
of reading activations at a specific layer in a neural network.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `X` | `ndarray` | Raw input. |
| `depth` | `int` | 0-indexed layer to stop at (inclusive). |

**Returns** Output of `layers_[depth].transform()`, with shape depending on the layer type.

**Raises** `IndexError` if `depth` is out of range.
**Raises** `RuntimeError` if called before `fit()`.

---

### BMCNetwork.predict

```python
BMCNetwork.predict(X) → ndarray, (M,)
```

Predict class labels using all layers.

**Returns** `ndarray, (M,), int32`

---

### BMCNetwork.predict_proba

```python
BMCNetwork.predict_proba(X) → ndarray, (M, C)
```

Return normalized probability distributions from the final layer.

**Returns** `ndarray, (M, C)` — each row sums to 1.0. Dtype follows input.

---

### BMCNetwork.predict_at

```python
BMCNetwork.predict_at(X, depth) → ndarray, (M,)
```

Return predictions using only layers `0..depth` (inclusive).

Enables depth-wise accuracy comparison without fitting multiple models.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `X` | `ndarray` | Raw input. |
| `depth` | `int` | Layer depth to stop at. |

**Returns** `ndarray, (M,), int32`

---

### BMCNetwork.__len__ / __getitem__

```python
len(net)     # number of layers
net[i]       # layer at position i
```

---

### BMCNetwork.get_state / set_state

```python
BMCNetwork.get_state() → dict
BMCNetwork.set_state(state, observation_spaces=None) → BMCNetwork
```

Serialize and restore network state.

`observation_spaces` — an `ObservationSpace` or list of them, one per layer that
requires one (`BMCLayer` layers). A single `ObservationSpace` is applied to layer 0;
the rest default to `None` (correct for `BMCCell` layers).

---

## bmc.composition.sequential — SequentialBMC

```python
class bmc.composition.sequential.SequentialBMC(
    bmc_layer,
    n_clusters_l2,
    metric_l2=None,
    interpolate_l2=False,
    rng=None,
)
```

Convenience constructor for the common 2-layer topology `BMCLayer → BMCCell`.

`SequentialBMC` is a subclass of `BMCNetwork`. All `BMCNetwork` methods are
available. This class adds `layer1_`, `layer2_` properties and `predict_l1()`
for backward compatibility.

Equivalent to:
```python
BMCNetwork([bmc_layer, BMCCell(n_clusters=n_clusters_l2, metric=metric_l2, ...)])
```

**Parameters**

| Name | Type | Default | Description |
|---|---|---|---|
| `bmc_layer` | `BMCLayer` | required | First layer. If already fitted, it will not be re-fitted. |
| `n_clusters_l2` | `int` | required | Number of clusters for the second-layer `BMCCell`. |
| `metric_l2` | `SimilarityMetric` | `L2Metric()` | Metric for the second layer. |
| `interpolate_l2` | `bool` | `False` | Geodesic interpolation for the second layer. |
| `rng` | `np.random.Generator` | `default_rng(42)` | |

**Properties**

| Name | Type | Description |
|---|---|---|
| `layer1_` | `BMCLayer` | The first layer. |
| `layer2_` | `BMCCell` | The second layer. |

**Additional methods**

```python
SequentialBMC.predict_l1(X) → ndarray, (M,)
```
Return predictions from layer 1 only — useful as a baseline to measure the
contribution of layer 2.

---

## bmc.composition.gate — BMCGate

```python
class bmc.composition.gate.BMCGate(
    n_clusters,
    interpolator=None,
    rng=None,
)
```

Geodesic-position router based on topology distortion. A BMCCell in
regression mode on 1D arc-length space. Learns the distortion landscape
along a parent cell's geodesic and produces routing weights.

The gate is mode-agnostic: the caller computes distortion (entropy for
classification, variance for regression) and passes it in.

**Parameters**

| Name | Type | Default | Description |
|---|---|---|---|
| `n_clusters` | `int` | required | Number of Voronoi cells for the gate's own geodesic on arc-length space. Typically set to the parent's K. |
| `interpolator` | `Interpolator` or `None` | `None` | Interpolation strategy for smooth routing. `LinearInterpolator()` recommended. `None` → hard (nearest gate cell). |
| `rng` | `np.random.Generator` | `default_rng(42)` | |

**Attributes** *(set after fit)*

| Name | Type | Description |
|---|---|---|
| `cell_` | `BMCCell` | Underlying regression BMCCell on 1D arc-length. |
| `distortion_` | `ndarray, (K_parent,)` | Per-cell distortion of the parent. |
| `max_distortion_` | `float` | Maximum distortion (used for normalization). |
| `is_fitted_` | `bool` | |

### BMCGate.fit

```python
BMCGate.fit(X, parent_cell, distortion) → BMCGate
```

Fit the gate on training data, a fitted parent cell, and precomputed
per-cell distortion values. Distortion is normalized to `[0, 1]` and
used as regression targets at each point's arc-length position.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `X` | `ndarray, (M, D)` | Training input features. |
| `parent_cell` | `BMCCell` | Fitted parent cell whose geodesic the gate operates on. |
| `distortion` | `ndarray, (K,), float64` | Per-cell distortion values. Zero = no information loss. |

### BMCGate.route

```python
BMCGate.route(X, parent_cell) → ndarray, (M,), float64
```

Compute routing weights for query points. Returns `w ∈ [0, 1]` per query:
`w ≈ 0` → topology preserved → prefer L1;
`w ≈ 1` → topology distorted → prefer L2.

If the gate has an interpolator, the weight varies smoothly along the
geodesic. Otherwise, it snaps to the nearest gate cell's value.

### BMCGate.get_state / set_state

```python
BMCGate.get_state() → dict
BMCGate.set_state(state) → BMCGate
```

---

## bmc.composition.skip — BMCSkipConnection

```python
class bmc.composition.skip.BMCSkipConnection(
    layer1,
    n_clusters_l2=None,
    impurity_threshold=0.0,
    gate_interpolator="linear",
    metric=None,
    rng=None,
)
```

Adaptive skip connection for resolving representation collapse. When a
BMCCell has impure Voronoi cells, a skip connection adds a specialist
Layer 2 trained on the impure regions, gated by a `BMCGate` on
Layer 1's geodesic arc-length.

**Entropy-weighted seeding:** L2's kmeans++ seed selection is weighted by
each training point's L1 cell distortion (entropy for classification,
variance for regression). Seeds concentrate in high-distortion regions,
greedy-allocating L2's limited K budget where L1 suffers most.

Works for both classification and regression via separate fit methods,
following the same pattern as `BMCCell`.

**Parameters**

| Name | Type | Default | Description |
|---|---|---|---|
| `layer1` | `BMCCell` | required | Primary cell. Can be pre-fitted or unfitted. |
| `n_clusters_l2` | `int` or `None` | `None` | K for the specialist Layer 2. If `None`, defaults to the number of training points in impure cells. |
| `n_clusters_gate` | `int` | `32` | K for the gate BMCCell. The gate operates on 1D arc-length — small K (16–64) is sufficient. Clamped to `min(n_clusters_gate, M)` at fit time. |
| `impurity_threshold` | `float` | `0.0` | Cells with impurity above this threshold are marked impure. For classification, 0.0 means any entropy = impure. For regression, use a small positive value (e.g., 0.1) to filter low-variance cells. |
| `gate_interpolator` | `Interpolator`, `"linear"`, or `None` | `"linear"` | Interpolation strategy for the gate. `"linear"` (default) → `LinearInterpolator()` for smooth routing. `None` → hard gating (snap to nearest gate cell). Pass any `Interpolator` instance for custom behavior. |
| `metric` | `SimilarityMetric` | `L2Metric()` | Metric for Layer 2. |
| `rng` | `np.random.Generator` | `default_rng(42)` | |

**Attributes** *(set after fit)*

| Name | Type | Description |
|---|---|---|
| `layer1_` | `BMCCell` | Fitted primary cell (frozen). |
| `gate_` | `BMCGate` or `None` | Gate on L1's geodesic arc-length. `None` if no impure cells. |
| `layer2_` | `BMCCell` or `None` | Specialist cell (entropy-weighted seeds). `None` if no impure cells. |
| `has_skip_` | `bool` | Whether a skip connection was needed. |
| `n_impure_cells_` | `int` | Number of impure cells in L1. |
| `n_routed_points_` | `int` | Number of training points routed to L2. |
| `impure_mask_` | `ndarray, (K,), bool` | Which L1 cells are impure. |
| `cell_impurity_` | `ndarray, (K,), float64` | Per-cell impurity values. |

---

### BMCSkipConnection.fit

```python
BMCSkipConnection.fit(X, y, n_classes=None) → BMCSkipConnection
```

Fit for classification. Impurity is measured as Shannon entropy of each
cell's class distribution. L2's seeds are weighted by per-cell entropy
so that higher-entropy regions receive denser seed coverage.

---

### BMCSkipConnection.fit_regression

```python
BMCSkipConnection.fit_regression(X, y) → BMCSkipConnection
```

Fit for regression. Impurity is measured as within-cell variance of
targets, normalized by global variance. L2's seeds are weighted by
per-cell variance ratio.

---

### BMCSkipConnection.predict / predict_regression

```python
BMCSkipConnection.predict(X) → ndarray, (M,), int32
BMCSkipConnection.predict_regression(X) → ndarray, (M,)
```

Predict with skip connection routing. The gate blends L1 and L2
predictions based on the learned distortion landscape along the geodesic.

---

### BMCSkipConnection.predict_l1 / predict_l1_regression

```python
BMCSkipConnection.predict_l1(X) → ndarray, (M,), int32
BMCSkipConnection.predict_l1_regression(X) → ndarray, (M,)
```

Return L1-only predictions (baseline, no skip).

---

### BMCSkipConnection.transform / predict_proba

```python
BMCSkipConnection.transform(X) → ndarray, (M, C)
BMCSkipConnection.predict_proba(X) → ndarray, (M, C)
```

Return probability distributions (classification mode).
`predict_proba` normalizes rows to sum to 1.

---

### Diagnostics

```python
BMCSkipConnection.entropy_per_cell() → ndarray, (K,)     # classification
BMCSkipConnection.impurity_per_cell() → ndarray, (K,)     # both modes
BMCSkipConnection.residual_entropy() → float              # classification
BMCSkipConnection.residual_impurity() → float             # both modes
BMCSkipConnection.gate_weights(X) → ndarray, (M,)        # routing weights per query
```

---

### BMCSkipConnection.get_state / set_state

```python
BMCSkipConnection.get_state() → dict
BMCSkipConnection.set_state(state) → BMCSkipConnection
```

Serialize and restore skip connection state including L1, gate, and L2.

---

## bmc.core.rvq — Residual Vector Quantization

```python
bmc.core.rvq.fit_rvq(cells, X, convergence_threshold=0.1) → int
```

Train a sequence of BMCCells using Residual Vector Quantization. Each cell
is fitted on the current residual via `fit_regression` (one pass per output
dimension). Returns the number of cells actually fitted.

**Parameters**

| Name | Type | Default | Description |
|---|---|---|---|
| `cells` | `list of BMCCell` | required | Pre-constructed cells. Fitted in place. |
| `X` | `ndarray, (M, D)` | required | Training data to reconstruct. |
| `convergence_threshold` | `float` | `0.1` | Stop when residual MSE < this. |

**Returns** `int` — number of cells used (may be < `len(cells)`).

---

```python
bmc.core.rvq.reconstruct_rvq(cells, X) → ndarray, (M, D)
```

Compute accumulated reconstruction from RVQ-fitted cells.

---

```python
bmc.core.rvq.residual_mse_history(cells, X) → list of float
```

Residual MSE after each layer.

---

## bmc.composition.rvq — SequentialRVQ

```python
class bmc.composition.rvq.SequentialRVQ(
    max_layers=64,
    k_per_layer=2,
    convergence_threshold=0.1,
    n_clusters_head=None,
    interpolator=None,
    metric=None,
    rng=None,
)
```

Convenience wrapper for sequential BMC trained via RVQ. Creates cells internally,
trains them with `fit_rvq()`, and attaches a prediction head.

**Parameters**

| Name | Type | Default | Description |
|---|---|---|---|
| `max_layers` | `int` | `64` | Maximum RVQ layers. May stop earlier. |
| `k_per_layer` | `int` | `2` | K for each layer. K=2 minimum for meaningful interpolation. |
| `convergence_threshold` | `float` | `0.1` | Stop when residual MSE < this. |
| `n_clusters_head` | `int` or `None` | `None` | K for prediction head. Defaults to `n_layers_used * k_per_layer`. |
| `interpolator` | `Interpolator` | `LinearInterpolator` | Template for regression cells. |
| `metric` | `SimilarityMetric` | `L2Metric()` | |
| `rng` | `np.random.Generator` | `default_rng(42)` | |

**Methods**

```python
SequentialRVQ.fit(X, y, n_classes=None) → SequentialRVQ           # classification
SequentialRVQ.fit_regression(X, y) → SequentialRVQ                # regression
SequentialRVQ.predict(X) → ndarray, (M,), int32                  # class labels
SequentialRVQ.predict_regression(X) → ndarray, (M,)              # continuous values
SequentialRVQ.predict_proba(X) → ndarray, (M, C)                 # probabilities
SequentialRVQ.transform(X) → ndarray, (M, D)                     # reconstruction
SequentialRVQ.reconstruction_error(X) → ndarray, (M,)            # per-sample MSE
SequentialRVQ.residual_history(X) → list of float                 # MSE per layer
SequentialRVQ.total_cells() → int                                 # total BMCCells
SequentialRVQ.total_k() → int                                     # total Voronoi cells
```

**Attributes** *(after fit)*: `n_layers_used_`, `n_classes_`, `head_`, `cells_`, `is_fitted_`

---

## bmc.metrics — Similarity Metrics

All metrics are first-class objects passed wherever a distance function is needed.
This mirrors PyTorch's design for loss functions.

All metrics accept an optional `device` parameter to control whether pairwise
computation runs on CPU (NumPy) or GPU (PyTorch CUDA). See
[`bmc.backend`](#bmcbackend--device-backend) for details.

### SimilarityMetric (abstract base)

```python
class bmc.metrics.base.SimilarityMetric(device=None)
```

Abstract base class. All subclasses must implement `pairwise()`.

**Parameters**

| Name | Type | Default | Description |
|---|---|---|---|
| `device` | `str` or `None` | `None` | Compute device for pairwise operations. `"cpu"` (default), `"cuda"` (requires PyTorch + CUDA GPU), or `"auto"` (use GPU if available, else CPU). `None` is equivalent to `"cpu"`. |

**Abstract method**

```python
SimilarityMetric.pairwise(A, B) → ndarray, (M, K)
```

Compute pairwise distances between all rows of A and all rows of B.

| Parameter | Type | Description |
|---|---|---|
| `A` | `ndarray, (M, D)` | M feature vectors. |
| `B` | `ndarray, (K, D)` | K reference vectors (typically seeds). |

**Returns** `dists` where `dists[i, j]` = distance between `A[i]` and `B[j]`.

**Convenience method**

```python
metric(a, b) → float
```

Compute distance between two individual vectors. Wraps `pairwise`.

**Device control**

```python
metric.to(device) → SimilarityMetric
```

Switch the metric to a different device. Returns `self` for chaining.

```python
metric.device → str
```

Read-only property returning the resolved device string (`"cpu"` or `"cuda"`).

---

### L2Metric

```python
class bmc.metrics.l2.L2Metric(device=None)
```

Euclidean (L2) distance between flat feature vectors.

$$d(a, b) = \|a - b\|_2 = \sqrt{\sum_i (a_i - b_i)^2}$$

Computed using the expand trick `||a||² + ||b||² - 2 a·b`. On CPU uses float64
to avoid precision loss. On CUDA runs in float32 on GPU via PyTorch. Returns float32.

**Parameters**

| Name | Type | Default | Description |
|---|---|---|---|
| `device` | `str` or `None` | `None` | `"cpu"`, `"cuda"`, or `"auto"`. |

**When to use:** General purpose. Default for all BMC components.

**GPU example:**
```python
metric = L2Metric(device="cuda")
cell = BMCCell(n_clusters=256, metric=metric)
cell.fit(X_train, y_train)  # all pairwise distances computed on GPU
```

---

### FrobeniusMetric

```python
class bmc.metrics.frobenius.FrobeniusMetric(device=None)
```

Frobenius norm distance between matrix-shaped features.

$$d(A, B) = \|A - B\|_F = \sqrt{\sum_{i,j} (A_{ij} - B_{ij})^2}$$

Accepts inputs of shape `(M, H, W)` — no flattening before distance computation.
For flat `(M, D)` inputs, reduces to `L2Metric`.

**Parameters**

| Name | Type | Default | Description |
|---|---|---|---|
| `device` | `str` or `None` | `None` | `"cpu"`, `"cuda"`, or `"auto"`. |

**When to use:** When features are naturally matrix-shaped (e.g. image patches,
covariance matrices) and element-wise differences carry geometric meaning.

---

### CosineMetric

```python
class bmc.metrics.cosine.CosineMetric(eps=1e-8, device=None)
```

Cosine distance: `1 - cosine_similarity`. Range `[0, 2]`.

$$d(a, b) = 1 - \frac{a \cdot b}{\|a\| \cdot \|b\|}$$

**Parameters**

| Name | Type | Default | Description |
|---|---|---|---|
| `eps` | `float` | `1e-8` | Added to norms before division to prevent divide-by-zero. |
| `device` | `str` or `None` | `None` | `"cpu"`, `"cuda"`, or `"auto"`. |

**When to use:** When vector magnitude is not meaningful — only direction matters.
Common for text embeddings and normalized feature spaces.

---

### SpectralMetric

```python
class bmc.metrics.spectral.SpectralMetric(device=None)
```

Spectral norm distance: largest singular value of the difference matrix.

$$d(A, B) = \|A - B\|_2 = \sigma_{\max}(A - B)$$

Requires computing an SVD per pair. O(D³) per element. Slower than L2/Frobenius
for large D, but sensitive to the dominant direction of variation.

**Parameters**

| Name | Type | Default | Description |
|---|---|---|---|
| `device` | `str` or `None` | `None` | Accepted for API consistency. Spectral norm requires per-pair SVD and currently runs on CPU regardless. |

**When to use:** When the principal component of the patch difference is the
meaningful signal — e.g. detecting dominant spatial frequencies in image patches.

---

### NuclearMetric

```python
class bmc.metrics.nuclear.NuclearMetric(device=None)
```

Nuclear norm distance: sum of all singular values of the difference matrix.

$$d(A, B) = \|A - B\|_* = \sum_i \sigma_i(A - B)$$

The convex envelope of matrix rank. Captures the total spread across all
directions, not just the dominant one. Dual norm of the spectral norm.

**Parameters**

| Name | Type | Default | Description |
|---|---|---|---|
| `device` | `str` or `None` | `None` | Accepted for API consistency. Nuclear norm requires per-pair SVD and currently runs on CPU regardless. |

**When to use:** When the full singular value spectrum of patch differences
is meaningful — e.g. capturing variation across all spatial frequency bands.

---

## bmc.observation — Observation Spaces

An observation space defines how raw input data is decomposed into zone feature
vectors. It is the only component the user must supply when working with non-tabular
data. Implementing a custom observation space requires subclassing `ObservationSpace`
and implementing one method: `extract()`.

### ObservationSpace (abstract base)

```python
class bmc.observation.base.ObservationSpace
```

**Abstract properties**

| Property | Type | Description |
|---|---|---|
| `n_zones` | `int` | Number of zones (parallel channels) this space produces. |
| `feature_dim` | `int` | Dimensionality D of each zone's feature vector. |

**Abstract method**

```python
ObservationSpace.extract(inputs) → ndarray, (N_zones, M, D)
```

Extract zone feature vectors from a batch of inputs.

| Parameter | Type | Description |
|---|---|---|
| `inputs` | `ndarray` | Batch of M raw inputs. Shape depends on domain. |

**Returns** `zone_vectors` where `zone_vectors[z, i]` is the D-dimensional feature
vector for sample i in zone z.

**Convenience**

```python
obs(inputs)  # equivalent to obs.extract(inputs)
```

---

### ImageObservationSpace

```python
class bmc.observation.image.ImageObservationSpace(
    image_size=28,
    pool_size=10,
    grid_n=5,
    output_size=5,
    method='lanczos',
    n_channels=1,
)
```

Decomposes a 2D image into a grid of overlapping patches, each pooled
to a fixed-size feature vector. Supports single-channel `(M, H, W)` and
multi-channel `(M, H, W, C)` inputs.

Multi-channel mode is used when one layer's spatial output feeds the next —
e.g., a vote matrix reshaped to `(M, grid_n, grid_n, C)`.

Grid layout: `grid_n × grid_n` patches extracted at evenly-spaced positions.
Position i on one axis: `round(i × (image_size - pool_size) / (grid_n - 1))`.

**Parameters**

| Name | Type | Default | Description |
|---|---|---|---|
| `image_size` | `int` | `28` | Height and width of input images (assumed square). |
| `pool_size` | `int` | `10` | Side length of each extracted patch. Must be ≤ `image_size`. |
| `grid_n` | `int` | `5` | Number of grid positions per axis. Total zones = `grid_n²`. Must be ≥ 2. |
| `output_size` | `int` | `5` | Side length of each pooled feature vector. |
| `method` | `str` | `'lanczos'` | Pooling method for single-channel inputs. See table below. |
| `n_channels` | `int` | `1` | Number of channels. Set to C for multi-channel input `(M, H, W, C)`. When > 1, patches include all channels: D = `output_size² × n_channels`. |

**Pooling methods** (single-channel only)

| Method | Backend | Notes |
|---|---|---|
| `'lanczos'` | PIL | High quality, no GPU. Recommended default. |
| `'bilinear'` | PyTorch | Requires `torch`. |
| `'bicubic'` | PyTorch | Requires `torch`. |
| `'avg'` | PyTorch | Adaptive average pooling. Requires `torch`. |
| `'max'` | PyTorch | Adaptive max pooling. Requires `torch`. |
| `'area'` | OpenCV | Anti-aliased. Requires `opencv-python`. |
| `'gaussian'` | OpenCV | Gaussian blur then area resize. Requires `opencv-python`. |

For multi-channel inputs: when `output_size == pool_size`, no pooling is applied
(patches used directly). Otherwise, spatial block averaging is used per channel.

**Properties**

| Property | Value |
|---|---|
| `n_zones` | `grid_n²` |
| `feature_dim` | `output_size²` (single-channel) or `output_size² × n_channels` (multi-channel) |

**extract() input shape**

- Single-channel: `(M, H, W)` or `(M, 1, H, W)`
- Multi-channel: `(M, H, W, C)` where `C == n_channels`

**Example: two-layer spatial BMC**

```python
# L1: 28x28 images → 5x5 grid of 7x7 patches
obs1 = ImageObservationSpace(image_size=28, pool_size=7, grid_n=5, output_size=7)
l1 = BMCLayer(obs1, n_clusters=1024)
l1.fit(images, labels, n_classes=62)
votes = l1.transform(images)  # (M, 25, 62)

# Reshape to spatial image: (M, 5, 5, 62)
vote_img = votes.reshape(M, 5, 5, 62)

# L2: 3x3 grid of 3x3 patches on the vote image
obs2 = ImageObservationSpace(image_size=5, pool_size=3, grid_n=3,
                              output_size=3, n_channels=62)
l2 = BMCLayer(obs2, n_clusters=256)
l2.fit(vote_img, labels, n_classes=62)  # 9 zones × 558D
```

**Raises** `ValueError` if `method` is not in the supported list, or `grid_n < 2`,
or `pool_size > image_size`.

**Additional methods**

```python
obs.grid_starts() → list of int
```
Return the list of patch start positions along one axis. Length = `grid_n`.

---

### TabularObservationSpace

```python
class bmc.observation.tabular.TabularObservationSpace(
    n_features,
    n_zones=None,
    zone_indices=None,
)
```

Decomposes a flat feature vector into N_zones subsets, each becoming an
independent zone with its own LUT.

Exactly one of `n_zones` or `zone_indices` must be specified.

**Parameters**

| Name | Type | Default | Description |
|---|---|---|---|
| `n_features` | `int` | required | Total number of input features F. |
| `n_zones` | `int`, optional | `None` | Automatic equal partitioning into N_zones groups. |
| `zone_indices` | `list of list of int`, optional | `None` | Explicit feature index list per zone. Zones may overlap. |

**Properties**

| Property | Value |
|---|---|
| `n_zones` | N_zones |
| `feature_dim` | Size of the largest zone (smaller zones are zero-padded). |

**extract() input shape** — `(M, F)` or `(F,)` for a single sample.

**Raises** `ValueError` if neither or both of `n_zones`/`zone_indices` are given,
or if any index in `zone_indices` is out of range.

**Additional methods**

```python
obs.zone_feature_counts() → list of int
```
Return the number of features in each zone.

---

### TimeSeriesObservationSpace

```python
class bmc.observation.timeseries.TimeSeriesObservationSpace(
    series_length,
    window_size,
    stride=1,
    n_features=1,
)
```

Decomposes a time series into overlapping sliding windows. Each window
becomes a zone with a flattened feature vector.

**Parameters**

| Name | Type | Default | Description |
|---|---|---|---|
| `series_length` | `int` | required | Length T of each input time series. |
| `window_size` | `int` | required | Length of each window. Must be ≤ `series_length`. |
| `stride` | `int` | `1` | Step between consecutive window starts. `stride=1` gives maximum overlap. `stride=window_size` gives non-overlapping windows. |
| `n_features` | `int` | `1` | Features per time step. 1 for univariate, >1 for multivariate. |

**Properties**

| Property | Value |
|---|---|
| `n_zones` | `(series_length - window_size) // stride + 1` |
| `feature_dim` | `window_size × n_features` |

**extract() input shape** — `(M, T)` for univariate, `(M, T, F)` for multivariate.

**Raises** `ValueError` if `window_size > series_length` or `stride < 1`.

**Additional methods**

```python
obs.window_starts() → list of int
```
Return the list of window start positions.

---

### TextObservationSpace

```python
class bmc.observation.text.TextObservationSpace(
    seq_len,
    embed_dim,
    chunk_size,
    stride=None,
)
```

Decomposes a sequence of token embeddings into contiguous chunks.
Each chunk becomes a zone with a flattened feature vector.

Input must be pre-embedded — this class does not perform tokenization
or embedding computation.

**Parameters**

| Name | Type | Default | Description |
|---|---|---|---|
| `seq_len` | `int` | required | Number of tokens/positions in the embedding sequence. |
| `embed_dim` | `int` | required | Dimensionality of each token embedding. |
| `chunk_size` | `int` | required | Number of consecutive embeddings per zone. Must be ≤ `seq_len`. |
| `stride` | `int`, optional | `chunk_size` | Step between chunk starts. Defaults to `chunk_size` (non-overlapping). Set `stride < chunk_size` for overlapping zones. |

**Properties**

| Property | Value |
|---|---|
| `n_zones` | `(seq_len - chunk_size) // stride + 1` |
| `feature_dim` | `chunk_size × embed_dim` |

**extract() input shape** — `(M, seq_len, embed_dim)`.

**Raises** `ValueError` if `chunk_size > seq_len`, `stride < 1`, or input shape mismatches.

**Additional methods**

```python
obs.chunk_starts() → list of int
```
Return the list of chunk start positions.

---

## bmc.model — High-Level Models

### BMCModel (abstract base)

```python
class bmc.model.base.BMCModel
```

Abstract base class for high-level BMC models. `BMCClassifier` and `BMCRegressor`
inherit from this. Provides the sklearn-compatible interface contract.

**Abstract methods**

```python
BMCModel.fit(X, y) → BMCModel
BMCModel.predict(X) → ndarray, (M,)
BMCModel.predict_proba(X) → ndarray, (M, C)
```

All subclasses must implement these three methods with the signatures above.

---

### BMCClassifier

```python
class bmc.model.classifier.BMCClassifier(
    observation_space=None,
    n_clusters=128,
    n_clusters_l2="auto",
    metric=None,
    interpolate=False,
    rng=None,
)
```

sklearn-compatible classifier wrapping a `BMCLayer` and a prediction head
(`BMCCell` operating on the vote matrix). Provides the easy path — sensible
defaults, single `fit`/`predict` API.

**Parameters**

| Name | Type | Default | Description |
|---|---|---|---|
| `observation_space` | `ObservationSpace`, optional | `None` | How to decompose inputs. If `None` and input is 2D tabular, a `TabularObservationSpace` with `n_zones = min(n_features, 10)` is created automatically. Must be set explicitly for image, time series, or text inputs. |
| `n_clusters` | `int` | `128` | Voronoi cells per zone in layer 1. |
| `n_clusters_l2` | `int`, `"auto"`, or `None` | `"auto"` | Voronoi cells for the prediction head (layer 2). `"auto"` = same as `n_clusters`. Set to `None` to disable the prediction head and use raw zone voting. |
| `metric` | `SimilarityMetric` | `L2Metric()` | Applied to all layers. |
| `interpolate` | `bool` | `False` | Geodesic interpolation for the prediction head (layer 2). Layer 1 always uses interpolation internally to produce smooth representations for the prediction head. |
| `rng` | `np.random.Generator` | `default_rng(42)` | |

**Interpolation policy:** Layer 1 (zone cells) always uses `interpolate=True` regardless
of this parameter — smooth intermediate representations produce more precise patterns for
the prediction head. The `interpolate` parameter controls only the prediction head (layer 2).

**Methods**

```python
BMCClassifier.fit(X, y) → BMCClassifier
BMCClassifier.predict(X) → ndarray, (M,), int32
BMCClassifier.predict_proba(X) → ndarray, (M, C)
```

**Raises** `ValueError` if `observation_space` is `None` and input is not 2D.
**Raises** `RuntimeError` if `predict` or `predict_proba` called before `fit`.

---

### BMCRegressor

```python
class bmc.model.regressor.BMCRegressor(
    observation_space=None,
    n_clusters=128,
    n_clusters_l2="auto",
    metric=None,
    interpolate=True,
    interpolator=None,
    rng=None,
)
```

sklearn-compatible regressor. Uses the same `BMCCell` infrastructure as
`BMCClassifier` — the only differences are that cells use `fit_regression()`
(LUT stores mean target values, not probability distributions), interpolation is
enabled by default, and `predict()` returns continuous values instead of argmax.

For multi-zone inputs, each zone independently predicts via a `BMCCell` fitted
with `fit_regression()`. When `n_clusters_l2` is set (default: `"auto"`), a
prediction head `BMCCell` clusters the zone prediction matrix and learns which
combinations of zone predictions map to which target values.

**Parameters**

| Name | Type | Default | Description |
|---|---|---|---|
| `observation_space` | `ObservationSpace`, optional | `None` | How to decompose inputs. Auto-created for 2D tabular with `n_zones = min(n_features, 10)`. |
| `n_clusters` | `int` | `128` | Voronoi cells per zone. |
| `n_clusters_l2` | `int`, `"auto"`, or `None` | `"auto"` | Clusters for the prediction head. `"auto"` = same as `n_clusters`. Set to `None` to disable and average zone predictions directly. |
| `metric` | `SimilarityMetric` | `L2Metric()` | Distance metric for all cells. |
| `interpolate` | `bool` | `True` | Geodesic interpolation for the prediction head (layer 2). If `True` and no `interpolator` is set, uses `LinearInterpolator`. |
| `interpolator` | `Interpolator` or `RDDAInterpolator` | `None` | Explicit interpolation strategy for the prediction head. Overrides `interpolate` flag. |
| `rng` | `np.random.Generator` | `default_rng(42)` | |

**Interpolation policy:** Zone cells (layer 1) always use `LinearInterpolator`
internally — smooth continuous predictions feed the prediction head more precise
patterns. The `interpolate` and `interpolator` parameters control only the prediction
head (layer 2).

**Methods**

```python
BMCRegressor.fit(X, y) → BMCRegressor
    # y : ndarray, (M,), continuous float
BMCRegressor.predict(X) → ndarray, (M,)
    # returns continuous predicted values
BMCRegressor.predict_proba(X) → ndarray, (M, N_zones)
    # returns per-zone predictions as columns (for inspection)
```

**Raises** `ValueError` if `observation_space` is `None` and input is not 2D.
**Raises** `RuntimeError` if `predict` or `predict_proba` called before `fit`.

---

## bmc.model.persistence — Save / Load

```python
from bmc.model.persistence import save, load
```

Models are serialized to `.bmc` files — ZIP archives containing `metadata.json`
(hyperparameters and non-array state) and `arrays/*.npy` (numpy arrays). No pickle
is used. Files are version-stamped and safe to inspect with any ZIP tool.

**Note:** The `ObservationSpace` is not serialized. You must supply it when loading
any model that uses one (`BMCLayer`, `BMCNetwork`, `BMCClassifier`, `BMCRegressor`).

---

### save

```python
bmc.model.persistence.save(model, path)
```

Save a fitted model to a `.bmc` file.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `model` | fitted BMC model | Any fitted `BMCCell`, `BMCLayer`, `BMCNetwork`, `SequentialBMC`, `BMCClassifier`, or `BMCRegressor`. |
| `path` | `str` or `Path` | Destination path. The `.bmc` extension is conventional but not required. |

**Raises** `RuntimeError` if the model is not fitted.

---

### load

```python
bmc.model.persistence.load(path, observation_space=None)
```

Load a fitted model from a `.bmc` file.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `path` | `str` or `Path` | Path to the `.bmc` file. |
| `observation_space` | `ObservationSpace` or `list of ObservationSpace`, optional | Required when loading models containing a `BMCLayer`. For `BMCNetwork`, pass a list with one entry per layer (use `None` for `BMCCell` layers). |

**Returns** A fitted model of the same class that was saved.

**Raises** `ValueError` if a `BMCLayer` is present and `observation_space` is not provided.

---

## bmc.core — Low-Level Functions

These functions are the engine underlying `BMCCell`. They are exposed publicly for
researchers who want to build custom pipelines that bypass the object-oriented API.

### bmc.core.seeding

```python
from bmc.core.seeding import kmeanspp_init, assign_clusters
```

---

#### kmeanspp_init

```python
bmc.core.seeding.kmeanspp_init(vectors, k, rng, metric=None, max_iter=50, tol=1e-4, lloyd_batch=None, batch_k=None, sample_weights=None) → ndarray, (K, D)
```

Initialize K seeds via kmeans++ then refine with Lloyd's algorithm.

**Phase 1 (kmeans++ init):** Seeds are chosen sequentially with probability
proportional to squared distance from the nearest already-chosen seed,
optionally scaled by `sample_weights`.

**Phase 2 (Lloyd's iteration):** Seeds are iteratively refined by assigning
a sample of vectors to the nearest seed, recomputing each seed as its cluster mean,
and repeating until convergence (max seed shift < `tol`) or `max_iter` is reached.
When `lloyd_batch < M`, each iteration uses a random subset (mini-batch Lloyd's).
A final full-dataset pass ensures the returned seeds are exact cluster means.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `vectors` | `ndarray, (M, D)` | Data points to seed from. Dtype policy applies. |
| `k` | `int` | Number of seeds. Must be ≤ M. |
| `rng` | `np.random.Generator` | Random generator. |
| `metric` | `SimilarityMetric`, optional | Defaults to `L2Metric()`. |
| `max_iter` | `int` | Maximum Lloyd's iterations. Default 50. |
| `tol` | `float` | Convergence tolerance on max seed shift. Default 1e-4. |
| `lloyd_batch` | `int` or `None` | Samples per Lloyd iteration. `None` → `min(M, 10*K)` — at least 10 samples per cluster for stable centroid estimates. Full-batch when M ≤ 10×K. |
| `batch_k` | `int` or `None` | Seed-batch size for pairwise distance computation in `assign_clusters`. Forwarded to every `assign_clusters` call (Lloyd iterations and final full-dataset pass). `None` → auto-compute from available memory inside `assign_clusters`. |
| `sample_weights` | `ndarray, (M,)`, optional | Per-point importance weights for seed selection. The sampling probability becomes `weight[i] * dist[i]**2 / sum(weight * dist**2)`. Biases seed placement toward high-weight regions (e.g., high-entropy cells from a previous layer). `None` → uniform weights (standard kmeans++). |

**Returns** `ndarray, (K, D)` — converged seed vectors (cluster centroids).

**Raises** `ValueError` if `k > M` or `k <= 0`.

**Complexity** O(K × M) for init, O(max_iter × K × min(M, 10K)) for refinement. Typically converges in 5–15 iterations.

---

#### assign_clusters

```python
bmc.core.seeding.assign_clusters(vectors, seeds, metric, batch_k=None) → ndarray, (M,), int32
```

Assign each vector to its nearest seed (Voronoi assignment).

Processes seeds in batches of `batch_k`, computing pairwise distances for each
batch and keeping a running minimum. When `batch_k >= K` the full (M, K) pairwise
matrix is computed in one vectorized call. When `batch_k == 1` this degenerates to
the original streaming behaviour (O(M) workspace).

**Parameters**

| Name | Type | Description |
|---|---|---|
| `vectors` | `ndarray, (M, D)` | Vectors to assign. |
| `seeds` | `ndarray, (K, D)` | Seed vectors. |
| `metric` | `SimilarityMetric` | Distance metric. |
| `batch_k` | `int` or `None` | Seeds per batch. `None` auto-computes from 75% of available memory (system RAM for CPU metrics, GPU VRAM for CUDA metrics). User overrides are clamped to `[1, K]` and warned if exceeding 99% of available memory. |

**Returns** `assignments` where `assignments[i]` is the index of the nearest seed to `vectors[i]`.

---

### bmc.core.geodesic

```python
from bmc.core.geodesic import tsp_nearest_neighbor, two_opt, build_tsp_path, tour_length
```

---

#### tsp_nearest_neighbor

```python
bmc.core.geodesic.tsp_nearest_neighbor(dists_matrix, start=0) → list of int
```

Construct an open TSP tour using the nearest-neighbor greedy heuristic.

Starts at `start`, always moves to the nearest unvisited node. Returns a
Hamiltonian path visiting all K nodes exactly once.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `dists_matrix` | `ndarray, (K, K)` | Pairwise distance matrix. |
| `start` | `int` | Starting node index. Default 0. |

**Returns** `list of int`, length K — ordered node indices.

**Complexity** O(K²). Returns an open tour (does not return to `start`).

---

#### two_opt

```python
bmc.core.geodesic.two_opt(path, dists_matrix) → list of int
```

Improve a TSP path using 2-opt local search.

Iteratively reverses sub-paths when doing so reduces total tour length.
Terminates when no improving swap exists (local optimum).

**Parameters**

| Name | Type | Description |
|---|---|---|
| `path` | `list of int` | Initial tour. |
| `dists_matrix` | `ndarray, (K, K)` | Pairwise distance matrix. |

**Returns** `list of int` — improved tour with length ≤ input tour length.

Does not modify the input list.

**Complexity** O(K²) per pass. Typically 2–5 passes for K ≤ 256.

---

#### build_tsp_path

```python
bmc.core.geodesic.build_tsp_path(dists_matrix, start=0) → ndarray, (K,), int64
```

Build an optimized TSP path: nearest-neighbor initialization + 2-opt improvement.

This is the primary entry point for geodesic ordering. The returned path defines
the binary address space — seed at position i has binary address i.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `dists_matrix` | `ndarray, (K, K)` | Pairwise distances between K seeds. |
| `start` | `int` | Starting node. Default 0. |

**Returns** `ndarray, (K,)` — TSP-ordered seed indices.

---

#### tour_length

```python
bmc.core.geodesic.tour_length(path, dists_matrix) → float
```

Compute the total path length (sum of consecutive edge distances).

**Parameters**

| Name | Type | Description |
|---|---|---|
| `path` | `array-like, length K` | Ordered node indices. |
| `dists_matrix` | `ndarray, (K, K)` | Pairwise distances. |

**Returns** `float` — total path length. Returns 0.0 for K ≤ 1.

---

### bmc.core.lut

```python
from bmc.core.lut import build_lut, build_lut_regression, reorder_lut, build_lut2, lut_purity, lut_dominant_class
```

---

#### build_lut

```python
bmc.core.lut.build_lut(assignments, labels, n_clusters, n_classes, dtype=None) → ndarray, (K, C)
```

Build a lookup table from cluster assignments and labels.

For each cluster k:
```
lut[k, c] = count(labels == c AND assignments == k) / count(assignments == k)
```

Empty clusters receive a uniform distribution `[1/C, ..., 1/C]`.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `assignments` | `ndarray, (M,), int` | Cluster index per sample. Values in `[0, n_clusters)`. |
| `labels` | `ndarray, (M,), int` | Class label per sample. Values in `[0, n_classes)`. |
| `n_clusters` | `int` | K — total number of clusters. |
| `n_classes` | `int` | C — total number of classes. |
| `dtype` | `dtype` or `None` | Output dtype. `None` → float32. Set by `BMCCell.fit()` from the input's dtype. |

**Returns** `ndarray, (K, C)` — each row is a valid probability distribution summing to 1.0. Dtype matches input.

---

#### build_lut_regression

```python
bmc.core.lut.build_lut_regression(assignments, targets, n_clusters, dtype=None) → ndarray, (K, 1)
```

Build a regression lookup table from cluster assignments and continuous targets.

For each cluster k:
```
lut[k, 0] = mean(targets[assignments == k])
```

Empty clusters receive the global mean. Rows are NOT normalized — they store
raw continuous values.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `assignments` | `ndarray, (M,), int` | Cluster index per sample. |
| `targets` | `ndarray, (M,), float` | Continuous target values. |
| `n_clusters` | `int` | K — total number of clusters. |
| `dtype` | `dtype` or `None` | Output dtype. `None` → float32. Set by `BMCCell.fit_regression()` from the input's dtype. |

**Returns** `ndarray, (K, 1)` — mean target value per cluster. Dtype matches input.

---

#### reorder_lut

```python
bmc.core.lut.reorder_lut(lut, tsp_path) → ndarray, (K, C)
```

Reorder LUT rows by TSP geodesic path.

After reordering, `lut[i]` corresponds to the i-th seed on the geodesic.
Applying the inverse permutation recovers the original ordering.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `lut` | `ndarray, (K, C)` | Original LUT. |
| `tsp_path` | `ndarray, (K,), int` | TSP-ordered indices from `build_tsp_path`. |

**Returns** `ndarray, (K, C)` — reordered LUT.

---

#### build_lut2

```python
bmc.core.lut.build_lut2(vote_matrix, labels, n_clusters, metric=None, rng=None) → tuple
```

Build a second-layer LUT from the vote matrix. Clusters the flattened vote patterns
(from all parallel zones) and builds a LUT mapping voting-pattern clusters to class
distributions. Applies the full pipeline: seeding → assignment → LUT → TSP reorder.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `vote_matrix` | `ndarray, (M, N_zones * C)` or `(M, N_zones, C)` | Flattened or unflattened vote matrix. Automatically flattened if 3D. |
| `labels` | `ndarray, (M,), int` | Class labels. |
| `n_clusters` | `int` | Number of second-layer clusters K2. |
| `metric` | `SimilarityMetric`, optional | Defaults to `L2Metric()`. |
| `rng` | `np.random.Generator`, optional | Defaults to `default_rng(42)`. |

**Returns** `(seeds2, lut2, tsp_path2)` — tuple of:

| Name | Type | Description |
|---|---|---|
| `seeds2` | `ndarray, (K2, D2)` | Second-layer cluster seeds in TSP order. |
| `lut2` | `ndarray, (K2, C)` | Second-layer lookup table in TSP order. Each row sums to 1.0. |
| `tsp_path2` | `ndarray, (K2,)` | TSP ordering of second-layer seeds. |

---

#### lut_purity

```python
bmc.core.lut.lut_purity(lut) → ndarray, (K,)
```

Compute the purity of each LUT row. Purity = `max(lut[k])`.

Values in `[1/C, 1.0]`. Higher purity means the cell is more class-pure.

---

#### lut_dominant_class

```python
bmc.core.lut.lut_dominant_class(lut) → ndarray, (K,), int32
```

Return `argmax(lut[k])` for each row — the most probable class per cell.

---

### bmc.core.inference

```python
from bmc.core.inference import (
    lookup_single_zone, accumulate_zones, predict,
    forward, normalize_distributions
)
```

---

#### lookup_single_zone

```python
bmc.core.inference.lookup_single_zone(vectors, seeds, lut, metric, batch_k=None) → ndarray, (M, C)
```

Hard LUT lookup for one zone: find nearest seed per vector, return its LUT row.
Delegates to `assign_clusters` internally and inherits its batched memory management.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `vectors` | `ndarray, (M, D)` | Input vectors. |
| `seeds` | `ndarray, (K, D)` | Seed vectors (TSP-ordered). |
| `lut` | `ndarray, (K, C)` | Lookup table (TSP-ordered). |
| `metric` | `SimilarityMetric` | Distance metric. |
| `batch_k` | `int` or `None` | Seed-batch size. See `assign_clusters`. |

**Returns** `ndarray, (M, C)` — `result[i] = lut[argmin_k d(vectors[i], seeds[k])]`.

---

#### accumulate_zones

```python
bmc.core.inference.accumulate_zones(zone_distributions) → ndarray, (M, C)
```

Sum probability distributions across parallel zones (soft voting).

**Parameters**

| Name | Type | Description |
|---|---|---|
| `zone_distributions` | `ndarray, (N_zones, M, C)` | Per-zone distributions. |

**Returns** `ndarray, (M, C)` — unnormalized sum. Use `argmax` directly or `normalize_distributions` first.

---

#### predict

```python
bmc.core.inference.predict(accumulated) → ndarray, (M,), int32
```

Convert accumulated distributions to hard predictions via `argmax`.

---

#### forward

```python
bmc.core.inference.forward(
    zone_vectors, zone_seeds, zone_luts, metric, return_distributions=False
) → ndarray, (M,) or tuple
```

Full forward pass through a parallel BMC layer. Runs `lookup_single_zone` per zone,
accumulates via soft voting, and returns argmax predictions.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `zone_vectors` | `ndarray, (N_zones, M, D)` | Input vectors per zone. |
| `zone_seeds` | `ndarray, (N_zones, K, D)` | Seeds per zone. |
| `zone_luts` | `ndarray, (N_zones, K, C)` | LUT per zone. |
| `metric` | `SimilarityMetric` | Distance metric. |
| `return_distributions` | `bool` | If `True`, return `(predictions, accumulated)` tuple. |

**Returns** `ndarray, (M,), int32` — hard predictions. If `return_distributions=True`,
returns `(predictions, distributions)` where `distributions` is `ndarray, (M, C)`. Dtype follows input.

---

#### normalize_distributions

```python
bmc.core.inference.normalize_distributions(distributions) → ndarray, (M, C)
```

Row-normalize accumulated distributions to valid probability vectors.

Zero-sum rows (which should not occur in practice) are handled gracefully
by substituting 1.0 as the denominator.

---

### bmc.core.interpolators — Interpolation Strategies

```python
from bmc.core.interpolators import (
    LinearInterpolator,
    SplineInterpolator,
    KernelInterpolator,
    RDDAInterpolator,
)
```

Four interpolation strategies are available. The first three (`Linear`, `Spline`,
`Kernel`) operate in arc-length space along the geodesic and share the `Interpolator`
base class. `RDDAInterpolator` operates in the original feature space and has a
distinct interface.

---

#### Interpolator (abstract base)

```python
class bmc.core.interpolators.Interpolator
```

**Abstract methods**

```python
Interpolator.fit(t_values, y_values) → Interpolator
    # t_values : ndarray, (K,) — arc-length positions, normalized to [0, 1]
    # y_values : ndarray, (K,) or (K, D) — stored values at each prototype

Interpolator.evaluate(t) → ndarray
    # t : ndarray, (M,) — arc-length positions to evaluate at
    # Returns: (M,) or (M, D) — interpolated values

Interpolator.get_state() → dict
Interpolator.set_state(state) → Interpolator
```

---

#### LinearInterpolator

```python
class bmc.core.interpolators.LinearInterpolator()
```

Piecewise-linear interpolation between adjacent prototypes. For each query position,
finds the bracketing prototypes and linearly blends their values.

Fast, O(1) per query. No fit-time computation beyond storing data. Default for
`BMCRegressor`.

---

#### SplineInterpolator

```python
class bmc.core.interpolators.SplineInterpolator(bc_type='not-a-knot')
```

Cubic spline through all prototypes. C2 continuous. Fits `scipy.interpolate.CubicSpline`
at training time. Default for `BMCCell` when `interpolate=True`.

**Parameters**

| Name | Type | Default | Description |
|---|---|---|---|
| `bc_type` | `str` | `'not-a-knot'` | Boundary condition. Options: `'not-a-knot'`, `'clamped'`, `'natural'`, `'periodic'`. |

Falls back to `LinearInterpolator` when fewer than 4 prototypes or degenerate arc-lengths.

---

#### KernelInterpolator

```python
class bmc.core.interpolators.KernelInterpolator(method='idw', power=2.0)
```

Inverse distance weighting or radial basis function interpolation.

**Parameters**

| Name | Type | Default | Description |
|---|---|---|---|
| `method` | `str` | `'idw'` | `'idw'` (inverse distance weighting) or `'rbf'` (radial basis function via `scipy.interpolate.RBFInterpolator`). |
| `power` | `float` | `2.0` | Exponent for IDW weighting. Only used for `method='idw'`. |

**Raises** `ValueError` if `method` is not `'idw'` or `'rbf'`.

---

#### RDDAInterpolator

```python
class bmc.core.interpolators.RDDAInterpolator(
    sensitivity='vc',
    k_neighbors=None,
    rectify=True,
)
```

Rectified Discriminatory Delta-Adjust interpolation in feature space. Unlike the
arc-length interpolators, RDDA operates directly in the original feature space:

$$\hat{f}(x) = \sum_i w_i \cdot [f(p_i) + \nabla f(p_i) \cdot (x - p_i)]$$

where $\nabla f(p_i)$ is estimated from neighbor finite differences in feature space,
and $w_i$ are sensitivity-weighted distance kernels.

**This class does not inherit from `Interpolator`** because it has a different
interface: `fit(seeds, values, metric)` and `evaluate(queries, seeds, metric)` take
feature-space arrays, not arc-length scalars. `BMCCell` and `BMCRegressor` detect
this via `isinstance` and route raw vectors accordingly.

**Parameters**

| Name | Type | Default | Description |
|---|---|---|---|
| `sensitivity` | `str` | `'vc'` | Sensitivity estimation method: `'fixed'` (constant), `'sa'` (local variance), `'vc'` (gradient magnitude — best overall), `'hod'` (curvature). |
| `k_neighbors` | `int` or `None` | `None` | Nearest prototypes per query in feature space. `None` = all K (global). **Recommend k=5 for best results.** |
| `rectify` | `bool` | `True` | Clamp adjustments so they don't go below the base prototype value. |

**Methods**

```python
RDDAInterpolator.fit(seeds, values, metric) → RDDAInterpolator
    # seeds  : ndarray, (K, D) — seed positions in feature space
    # values : ndarray, (K,) or (K, C) — stored values
    # metric : SimilarityMetric

RDDAInterpolator.evaluate(queries, seeds, metric) → ndarray
    # queries : ndarray, (M, D) — raw query vectors
    # seeds   : ndarray, (K, D) — same seeds as fit
    # metric  : SimilarityMetric
    # Returns : (M,) or (M, C)

RDDAInterpolator.get_state() → dict
RDDAInterpolator.set_state(state) → RDDAInterpolator
```

**Raises** `ValueError` if `sensitivity` is not one of `'fixed'`, `'sa'`, `'vc'`, `'hod'`.

---

### bmc.core.interpolation

```python
from bmc.core.interpolation import (
    project_onto_segment,
    project_onto_segment_batch,
    interpolate_distributions,
    interpolate_distributions_batch,
    lookup_interpolated_single_zone,
    forward_interpolated,
)
```

---

#### project_onto_segment

```python
bmc.core.interpolation.project_onto_segment(query, p0, p1) → float
```

Project a query point onto the line segment `[p0, p1]`.

Returns scalar `t ∈ [0, 1]` such that the projection of `query` onto the
segment is `p0 + t × (p1 - p0)`. Values outside `[0, 1]` are clamped.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `query` | `ndarray, (D,)` | Query point. |
| `p0` | `ndarray, (D,)` | Segment start. |
| `p1` | `ndarray, (D,)` | Segment end. |

**Returns** `float` in `[0, 1]`. Returns 0.0 for degenerate segments (`p0 == p1`).

---

#### project_onto_segment_batch

```python
bmc.core.interpolation.project_onto_segment_batch(queries, p0, p1) → ndarray, (M,)
```

Vectorized version of `project_onto_segment` for a batch of queries.

---

#### interpolate_distributions

```python
bmc.core.interpolation.interpolate_distributions(lut_i, lut_j, t) → ndarray, (C,)
```

Linearly interpolate between two LUT rows.

$$f = (1 - t) \cdot \text{lut}_i + t \cdot \text{lut}_j$$

If both inputs are valid probability distributions (non-negative, sum to 1),
the output is guaranteed to be a valid probability distribution.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `lut_i` | `ndarray, (C,)` | Distribution at `t=0`. |
| `lut_j` | `ndarray, (C,)` | Distribution at `t=1`. |
| `t` | `float` | Interpolation parameter. Clamped to `[0, 1]`. |

---

#### interpolate_distributions_batch

```python
bmc.core.interpolation.interpolate_distributions_batch(lut_i, lut_j, t) → ndarray, (M, C)
```

Batch linear interpolation for M queries.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `lut_i` | `ndarray, (C,)` | Distribution at `t=0`. |
| `lut_j` | `ndarray, (C,)` | Distribution at `t=1`. |
| `t` | `ndarray, (M,)` | Per-query interpolation parameters. Clamped to `[0, 1]`. |

---

#### lookup_interpolated_single_zone

```python
bmc.core.interpolation.lookup_interpolated_single_zone(
    vectors, seeds, lut, tsp_path, metric
) → ndarray, (M, C)
```

Interpolated LUT lookup for a single zone.

For each query vector:
1. Find nearest seed `p_i` by Voronoi assignment
2. Select the closer TSP neighbor `j*` from `{(i-1) % K, (i+1) % K}` — symmetric selection ensures queries on either side of `p_i` blend in the correct direction
3. Project query onto segment `[p_i, p_{j*}]` to obtain `t` — projection-based blending measures position along the geodesic segment, preventing over-blending in high dimensions
4. Return `(1 - t) × lut[i] + t × lut[j*]`

**Note on interpolation and feature scaling:** Interpolation assumes the distance
metric is geometrically meaningful. With pathologically unscaled features (e.g. one
feature in `[0, 1000]` and another in `[0, 0.001]`), the geodesic and projection
are distorted. In such cases, either scale features before fitting or use hard
lookup (`interpolate=False`).

**Parameters**

| Name | Type | Description |
|---|---|---|
| `vectors` | `ndarray, (M, D)` | Input vectors. |
| `seeds` | `ndarray, (K, D)` | Seeds in TSP order. |
| `lut` | `ndarray, (K, C)` | LUT in TSP order. |
| `tsp_path` | `ndarray, (K,)` | TSP path (used to derive neighbor indices). |
| `metric` | `SimilarityMetric` | Distance metric. |

**Returns** `ndarray, (M, C)` — interpolated distributions. Each row sums to 1.0. Dtype follows input.

---

#### forward_interpolated

```python
bmc.core.interpolation.forward_interpolated(
    zone_vectors, zone_seeds, zone_luts, zone_tsp_paths, metric,
    return_distributions=False,
) → ndarray, (M,) or tuple
```

Full interpolated forward pass through a parallel BMC layer. Runs
`lookup_interpolated_single_zone` per zone, accumulates via soft voting,
and returns argmax predictions. The interpolated counterpart of
`bmc.core.inference.forward`.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `zone_vectors` | `ndarray, (N_zones, M, D)` | Input vectors per zone. |
| `zone_seeds` | `ndarray, (N_zones, K, D)` | Seeds in TSP order per zone. |
| `zone_luts` | `ndarray, (N_zones, K, C)` | LUTs in TSP order per zone. |
| `zone_tsp_paths` | `ndarray, (N_zones, K)` | TSP path per zone. |
| `metric` | `SimilarityMetric` | Distance metric. |
| `return_distributions` | `bool` | If `True`, return `(predictions, accumulated)` tuple. |

**Returns** `ndarray, (M,), int32` — hard predictions. If `return_distributions=True`,
returns `(predictions, distributions)` where `distributions` is `ndarray, (M, C)`. Dtype follows input.

---

## Implementing a Custom ObservationSpace

To use BMC on a domain not covered by the built-in spaces, subclass
`ObservationSpace` and implement `n_zones`, `feature_dim`, and `extract()`:

```python
import numpy as np
from bmc.observation.base import ObservationSpace

class GraphObservationSpace(ObservationSpace):
    """
    Example: decompose a graph into ego-network zones.
    Each zone is the feature vector of a node's local neighborhood.
    """
    def __init__(self, node_indices, neighbor_depth=1):
        self.node_indices = node_indices  # list of central nodes per zone
        self.neighbor_depth = neighbor_depth
        self._feature_dim = ...  # set based on your aggregation

    @property
    def n_zones(self) -> int:
        return len(self.node_indices)

    @property
    def feature_dim(self) -> int:
        return self._feature_dim

    def extract(self, inputs: np.ndarray) -> np.ndarray:
        # inputs: (M, ...) — your graph representation
        # returns: (N_zones, M, D)
        result = np.zeros((self.n_zones, len(inputs), self.feature_dim), dtype=np.float32)
        for z, node_idx in enumerate(self.node_indices):
            result[z] = self._aggregate_neighborhood(inputs, node_idx)
        return result

    def _aggregate_neighborhood(self, inputs, node_idx):
        ...  # domain-specific logic
```

Once implemented, it plugs directly into `BMCLayer` and `BMCClassifier` without
any changes to the rest of the framework:

```python
obs = GraphObservationSpace(node_indices=[0, 5, 12, 7])
layer = BMCLayer(obs, n_clusters=256)
layer.fit(graphs, labels)
```

---

## bmc.backend — Device Backend

```python
from bmc.backend import get_backend, resolve_device
```

Provides a thin abstraction over NumPy (CPU) and PyTorch (GPU) for the operations
that dominate training time: pairwise distances, argmin, and matrix multiply.

The backend is **not** a full array library — it only wraps the specific operations
that are hot in the BMC training loop. Everything else stays in NumPy.

### Installation

GPU support requires PyTorch with CUDA:

```bash
pip install bmc[gpu]          # installs torch>=2.0
# or
pip install biomimetic-framework[gpu]
```

---

### resolve_device

```python
bmc.backend.resolve_device(device) → str
```

Resolve a device string to an actual usable device.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `device` | `str` or `None` | `"cpu"`, `"cuda"`, `"auto"`, or `None`. `None` and `"auto"` both auto-detect: use CUDA if PyTorch is installed and a CUDA GPU is available, else CPU. |

**Returns** `"cpu"` or `"cuda"`.

**Raises** `RuntimeError` if `device="cuda"` but PyTorch is not installed or no CUDA GPU is available.
**Raises** `ValueError` if `device` is not a recognized string.

---

### get_backend

```python
bmc.backend.get_backend(device=None) → CPUBackend or CUDABackend
```

Return the backend class for the given device.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `device` | `str` or `None` | Passed to `resolve_device`. |

**Returns** `CPUBackend` or `CUDABackend` — a class with static methods for pairwise
distance computation.

---

### CPUBackend

NumPy-based backend. All methods are static and operate on NumPy arrays.

| Method | Signature | Description |
|---|---|---|
| `pairwise_l2` | `(A, B) → (M, K)` | L2 distance. CPU uses float64 internally, returns input dtype. |
| `pairwise_cosine` | `(A, B, eps=1e-8) → (M, K)` | Cosine distance. Returns input dtype. |
| `pairwise_frobenius` | `(A, B) → (M, K)` | Frobenius distance (flattens inputs). Returns input dtype. |
| `argmin` | `(dists, axis=1) → int32` | Argmin over distance matrix. |
| `bincount` | `(x, weights=None, minlength=0) → ndarray` | Weighted bincount. |

---

### CUDABackend

PyTorch CUDA backend. Same interface as `CPUBackend`. Accepts NumPy arrays, moves
them to GPU, runs the operation, and returns NumPy arrays. GPU tensors are freed
after each call.

For large `M × K` products (the pairwise distance matrix), this is dramatically
faster than CPU — matrix multiply is the canonical GPU workload.

| Method | Signature | Notes |
|---|---|---|
| `pairwise_l2` | `(A, B) → (M, K)` | Runs in float32 on GPU. Returns input dtype. |
| `pairwise_cosine` | `(A, B, eps=1e-8) → (M, K)` | Normalize + matmul on GPU. Returns input dtype. |
| `pairwise_frobenius` | `(A, B) → (M, K)` | Identical to `pairwise_l2` (flattened). Returns input dtype. |
| `argmin` | `(dists, axis=1) → int32` | `torch.argmin` on GPU. |
| `bincount` | `(x, weights=None, minlength=0) → ndarray` | Stays on CPU (not a GPU win). |

---

### Usage pattern

Pass a `device`-configured metric to any BMC component:

```python
from bmc import BMCCell, BMCLayer, BMCClassifier, L2Metric

# GPU-accelerated training
metric = L2Metric(device="cuda")
cell = BMCCell(n_clusters=256, metric=metric)
cell.fit(X_train, y_train)

# Auto-detect: use GPU if available, else CPU
clf = BMCClassifier(n_clusters=128, metric=L2Metric(device="auto"))
clf.fit(X_train, y_train)

# Switch an existing metric
metric = L2Metric()
metric.to("cuda")
```

All pairwise distance calls in kmeans++ init, Lloyd's iteration, TSP distance matrix
construction, and inference are routed through the backend. No other code changes
are needed.

---

## bmc.fpga — FPGA Export

Export a trained `BMCCell` to FPGA-ready files (Verilog RTL, ROM hex, testbench).

```python
from bmc import BMCCell, L2Metric, export_for_fpga

cell = BMCCell(n_clusters=128, metric=L2Metric())
cell.fit(X_train, y_train, n_classes=10)

result = export_for_fpga(cell, "build/", X_test=X_test, y_test=y_test)
print(result["quantized_accuracy"])  # e.g. 0.97
```

### export_for_fpga

```python
export_for_fpga(
    cell,
    output_dir: str,
    X_test: np.ndarray = None,
    y_test: np.ndarray = None,
    data_min: float = None,
    data_max: float = None,
) -> dict
```

**Parameters**

| Parameter | Type | Description |
|---|---|---|
| `cell` | `BMCCell` | A fitted cell (must have `seeds_` and `lut_`). |
| `output_dir` | `str` | Target directory. Created if it does not exist. |
| `X_test` | `np.ndarray` | Optional test vectors, shape `(N, D)`. Enables testbench generation and quantized accuracy verification. |
| `y_test` | `np.ndarray` | Optional test labels, shape `(N,)`. Required if `X_test` is given. |
| `data_min` | `float` | Minimum value for 8-bit quantization. Inferred from seeds and test data if omitted. |
| `data_max` | `float` | Maximum value for 8-bit quantization. Inferred from seeds and test data if omitted. |

**Returns**

A `dict` with:

| Key | Type | Description |
|---|---|---|
| `K` | `int` | Number of seeds (clusters). |
| `D` | `int` | Feature dimensions. |
| `C` | `int` | Number of classes. |
| `model_bytes` | `int` | Total ROM size (`K*D + K*C`). |
| `seeds_bytes` | `int` | Seeds ROM size. |
| `lut_bytes` | `int` | LUT ROM size. |
| `quantization_scale` | `float` | Scale factor used for 8-bit quantization. |
| `quantization_offset` | `float` | Offset used for 8-bit quantization. |
| `quantized_accuracy` | `float` | Accuracy using integer-only arithmetic (present if test data given). |
| `test_correct` | `int` | Number of correct predictions (present if test data given). |
| `test_total` | `int` | Total test samples (present if test data given). |

---

## bmc.viz — Visualization

Architecture summaries, diagrams, data scatter plots, and composite visualizations.

```python
from bmc.viz import summary, diagram, scatter, architecture
```

Optional dependencies: `pip install bmc[viz]` (installs `rich`, `graphviz`).
`scatter` and `architecture` additionally require `umap-learn` and `matplotlib`.

---

### summary

```python
bmc.viz.summary(model, return_text=False)
```

Print a CLI architecture tree showing components, parameters, and memory footprint.
Uses `rich` if installed, falls back to plain text.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `model` | any BMC model | BMCCell, BMCLayer, BMCNetwork, BMCClassifier, or BMCRegressor. |
| `return_text` | `bool` | If `True`, return the summary as a string instead of printing. |

---

### diagram

```python
bmc.viz.diagram(model, output_path="bmc_architecture", fmt=None)
```

Generate a Graphviz architecture diagram. Black and white, LTR horizontal layout.
Shows components, parameters, memory footprint. No fitted-state information
(purity, status). Format inferred from file extension.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `model` | any BMC model | |
| `output_path` | `str` | Output path. Defaults to `bmc_architecture.svg`. |
| `fmt` | `str`, optional | `"svg"`, `"pdf"`, or `"png"`. Inferred from extension if omitted. |

**Returns** `str` — path to the generated file.

---

### scatter

```python
bmc.viz.scatter(model, X, y, output_path="bmc_scatter", fmt=None,
                class_names=None, layer=None, zone=None, method="umap")
```

2D scatter plot with Voronoi regions and geodesic path. Projects data and seeds
to 2D via UMAP or PCA, draws Voronoi polygons colored by dominant class, and
shows the TSP geodesic path as directed arrows.

- **BMCCell**: single plot
- **BMCLayer**: grid of all zones, or single zone with `zone=N`
- **BMCNetwork**: selects a layer, then grid or single zone
- **BMCClassifier**: unwraps to internal network
- **BMCRegressor**: grid of zone cells, or single zone

**Parameters**

| Name | Type | Description |
|---|---|---|
| `model` | any BMC model | |
| `X` | `ndarray` | Data to plot (training data recommended). |
| `y` | `ndarray, (N,), int` | Integer labels for coloring. |
| `output_path` | `str` | Output file path. |
| `class_names` | `list of str`, optional | Legend labels. |
| `layer` | `int`, optional | Layer index for BMCNetwork. |
| `zone` | `int`, optional | Zone index for BMCLayer. |
| `method` | `str` | `"umap"` or `"pca"`. Default: `"umap"`. |

**Returns** `str` — path to the generated file.

---

### architecture

```python
bmc.viz.architecture(model, X, y, output_path="bmc_architecture",
                     fmt=None, class_names=None, method="umap")
```

Composite visualization: architecture pipeline with actual scatter plots
embedded inside each cell box. Replaces text labels with Voronoi + geodesic
plots. Shows observation space label, zone cell scatter grid, and prediction
head scatter in a single LTR figure.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `model` | any fitted BMC model | |
| `X` | `ndarray` | Training data. |
| `y` | `ndarray, (N,), int` | Integer labels for coloring. |
| `output_path` | `str` | Output file path. |
| `class_names` | `list of str`, optional | Legend labels. |
| `method` | `str` | `"umap"` or `"pca"`. Default: `"umap"`. |

**Returns** `str` — path to the generated file.

---

### trace

```python
bmc.viz.trace(model, X, y, query, output_path="bmc_trace",
              fmt=None, class_names=None, method="umap")
```

Prediction trace: show which Voronoi cells and geodesic edges activate for
a single query point. For each cell the query passes through: seed i gets
yellow highlight (alpha=t), seed j gets yellow highlight (alpha=1-t), the
active geodesic edge is highlighted gold, and a purple dot marks the
parametric position. Recursive across zones and layers.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `model` | any fitted BMC model | |
| `X` | `ndarray` | Training data (for background scatter). |
| `y` | `ndarray, (N,), int` | Integer labels. |
| `query` | `ndarray, (D,)` | Single query point. |
| `output_path` | `str` | Output file path. |
| `class_names` | `list of str`, optional | Legend labels. |
| `method` | `str` | `"umap"` or `"pca"`. Default: `"umap"`. |

**Returns** `str` — path to the generated file.

---

### parametric

```python
bmc.viz.parametric(cell, output_path="bmc_parametric", fmt=None,
                   class_names=None, n_samples=500)
```

Geodesic parametric plot: curvature (Y-axis) vs arc-length (X-axis), with
color-gradient blending from the interpolator's class predictions. Shows
where the geodesic bends sharply (class boundaries) and where it runs
straight (cluster interiors). Seed positions marked as dots.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `cell` | `BMCCell` (fitted) | |
| `output_path` | `str` | Output file path. |
| `class_names` | `list of str`, optional | Legend labels (shown for ≤10 classes). |
| `n_samples` | `int` | Dense sampling resolution along geodesic. Default: 500. |

**Returns** `str` — path to the generated file.

---

### riemannian

```python
bmc.viz.riemannian(cell, output_path="bmc_riemannian", fmt=None,
                   class_names=None, grid_samples=50000, padding=0.1)
```

3D Riemannian manifold: the BMC's learned function surface in (D+C)-dimensional
joint input-output space. Grids the input space, predicts output at each point,
and renders the manifold as a colored mesh with the geodesic curve and seed
markers. For D+C ≤ 3, plotted directly. For D+C > 3, projected via UMAP.

**Parameters**

| Name | Type | Description |
|---|---|---|
| `cell` | `BMCCell` (fitted) | |
| `output_path` | `str` | Output file path. |
| `class_names` | `list of str`, optional | |
| `grid_samples` | `int` | Number of grid points to sample. Default: 50000. |
| `padding` | `float` | Fraction of input range to pad for extrapolation. Default: 0.1. |

**Returns** `str` — path to the generated file.

---

## Exceptions

| Exception | Raised by | Condition |
|---|---|---|
| `RuntimeError` | `BMCCell.transform`, `BMCCell.predict`, `BMCCell.purity`, `BMCCell.dominant_classes`, `BMCCell.get_state` | Called before `fit()`. |
| `RuntimeError` | `BMCLayer.transform`, `BMCLayer.predict`, `BMCLayer.predict_proba` | Called before `fit()`. |
| `RuntimeError` | `BMCNetwork.transform`, `BMCNetwork.predict`, `BMCNetwork.predict_proba`, `BMCNetwork.transform_at`, `BMCNetwork.predict_at` | Called before `fit()`. |
| `RuntimeError` | `BMCClassifier.predict`, `BMCClassifier.predict_proba` | Called before `fit()`. |
| `RuntimeError` | `BMCRegressor.predict`, `BMCRegressor.predict_proba` | Called before `fit()`. |
| `RuntimeError` | `save()` | Model is not fitted. |
| `ValueError` | `BMCNetwork.__init__` | `layers` is empty. |
| `ValueError` | `BMCNetwork.set_state` | A `BMCLayer` layer has no `ObservationSpace` supplied. |
| `ValueError` | `kmeanspp_init` | `k > M` or `k <= 0`. |
| `ValueError` | `ImageObservationSpace.__init__` | Unsupported `method`, `grid_n < 2`, or `pool_size > image_size`. |
| `ValueError` | `TabularObservationSpace.__init__` | Neither or both of `n_zones`/`zone_indices` given; invalid index in `zone_indices`. |
| `ValueError` | `TimeSeriesObservationSpace.__init__` | `window_size > series_length` or `stride < 1`. |
| `ValueError` | `TextObservationSpace.__init__` | `chunk_size > seq_len` or `stride < 1`. |
| `TypeError` | `BMCNetwork.__init__` | A layer is missing `fit`, `transform`, or `predict`. |
| `IndexError` | `BMCNetwork.transform_at`, `BMCNetwork.predict_at` | `depth` out of range. |
| `ValueError` | `KernelInterpolator.__init__` | `method` not `'idw'` or `'rbf'`. |
| `ValueError` | `RDDAInterpolator.__init__` | `sensitivity` not one of `'fixed'`, `'sa'`, `'vc'`, `'hod'`. |
| `RuntimeError` | `resolve_device` | `device="cuda"` but PyTorch is not installed or no CUDA GPU is available. |
| `ValueError` | `resolve_device` | Unrecognized `device` string. |
| `ImportError` | `ImageObservationSpace._pool_*` | Required backend (Pillow, PyTorch, OpenCV) not installed. |
| `ValueError` | `export_for_fpga` | Cell is not fitted (no `seeds_` attribute). |
| `RuntimeError` | `BMCCell.predict_regression` | Called before `fit_regression()`. |
| `ValueError` | `scatter`, `architecture` | Model is not fitted, or `method` not `"umap"` or `"pca"`. |
| `ImportError` | `scatter`, `architecture` | `umap-learn` not installed. |
| `ImportError` | `diagram` | `graphviz` not installed. |
| `ValueError` | `trace`, `parametric`, `riemannian` | Model/cell is not fitted. |
| `ImportError` | `riemannian` | `umap-learn` not installed (when D+C > 3). |
