biomimetic-framework 0.1.3

Biomimetic Lookup Cell — parametric geodesic lookup framework for machine learning

BMC Framework — API Reference

Version 0.1.2

---

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:

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

# Composition
from bmc import BMCLayer, BMCNetwork, SequentialBMC

# 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

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

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

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

---

bmc.core.cell — BMCCell

class bmc.core.cell.BMCCell(
    n_clusters,
    metric=None,
    interpolate=False,
    interpolator=None,
    rng=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.

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

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

Fit the cell on training vectors and labels.

Steps executed:

1. Initialize K seeds via kmeans++ from vectors
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), float32	Training feature vectors for this zone.
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.

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).

---

BMCCell.transform

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), float64	Per-sample probability distributions. Each row sums to 1.0.

Raises RuntimeError if called before fit().

---

BMCCell.predict

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.purity

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,), float64	Per-cell purity values in [1/C, 1.0].

Raises RuntimeError if called before fit().

---

BMCCell.dominant_classes

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

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

class bmc.composition.layer.BMCLayer(
    observation_space,
    n_clusters,
    metric=None,
    interpolate=False,
    rng=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.

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

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

Fit all zone cells on training data.

Zone fitting is parallelized across CPU cores via ProcessPoolExecutor when there are 4 or more zones. Each zone's BMCCell is fitted independently on its own zone vectors. 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

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), float64	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.

Raises RuntimeError if called before fit().

---

BMCLayer.predict

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

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

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

Returns

Name	Type	Description
probabilities	ndarray, (M, C), float64	Each row sums to 1.0.

---

bmc.composition.sequential — BMCNetwork

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

# 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

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

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

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

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

Predict class labels using all layers.

Returns ndarray, (M,), int32

---

BMCNetwork.predict_proba

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

Return normalized probability distributions from the final layer.

Returns ndarray, (M, C), float64 — each row sums to 1.0.

---

BMCNetwork.predict_at

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

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

---

BMCNetwork.get_state / set_state

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

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:

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

SequentialBMC.predict_l1(X) → ndarray, (M,)

Return predictions from layer 1 only — useful as a baseline to measure the contribution of layer 2.

---

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 for details.

SimilarityMetric (abstract base)

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

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

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

metric(a, b) → float

Compute distance between two individual vectors. Wraps pairwise.

Device control

metric.to(device) → SimilarityMetric

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

metric.device → str

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

---

L2Metric

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:

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

---

FrobeniusMetric

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

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

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

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)

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

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

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

---

ImageObservationSpace

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

Decomposes a 2D grayscale image into a grid of overlapping patches, each pooled to a fixed-size feature vector.

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 in pixels. 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. Feature dimension D = output_size².
method	str	'lanczos'	Pooling method. See table below.

Pooling methods

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.

Properties

Property	Value
n_zones	grid_n²
feature_dim	output_size²

extract() input shape — (M, H, W) or (M, 1, H, W), values in [0, 1].

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

Additional methods

obs.grid_starts() → list of int

Return the list of patch start positions along one axis. Length = grid_n.

---

TabularObservationSpace

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

obs.zone_feature_counts() → list of int

Return the number of features in each zone.

---

TimeSeriesObservationSpace

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

obs.window_starts() → list of int

Return the list of window start positions.

---

TextObservationSpace

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

obs.chunk_starts() → list of int

Return the list of chunk start positions.

---

bmc.model — High-Level Models

BMCModel (abstract base)

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

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

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

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

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

---

BMCRegressor

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. Stores the mean target value at each Voronoi cell and uses geodesic interpolation to predict continuous values. No bin discretization.

For multi-zone inputs, each zone independently predicts via its own regression cell. When n_clusters_l2 is set (default: "auto"), a prediction head clusters the zone prediction matrix and learns which combinations of zone predictions map to which target values — analogous to the classifier's LUT2.

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

BMCRegressor.fit(X, y) → BMCRegressor
    # y : ndarray, (M,), continuous float
BMCRegressor.predict(X) → ndarray, (M,), float64
    # returns continuous predicted values
BMCRegressor.predict_proba(X) → ndarray, (M, N_zones), float64
    # 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

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

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

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

from bmc.core.seeding import kmeanspp_init, assign_clusters

---

kmeanspp_init

bmc.core.seeding.kmeanspp_init(vectors, k, rng, metric=None, max_iter=50, tol=1e-4) → 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.

Phase 2 (Lloyd's iteration): Seeds are iteratively refined by assigning all 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.

Parameters

Name	Type	Description
vectors	ndarray, (M, D), float32	Data points to seed from.
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.

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 × M) for refinement. Typically converges in 5–15 iterations.

---

assign_clusters

bmc.core.seeding.assign_clusters(vectors, seeds, metric) → ndarray, (M,), int32

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

Parameters

Name	Type	Description
vectors	ndarray, (M, D)	Vectors to assign.
seeds	ndarray, (K, D)	Seed vectors.
metric	SimilarityMetric	Distance metric.

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

---

bmc.core.geodesic

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

---

tsp_nearest_neighbor

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

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

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

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

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

---

build_lut

bmc.core.lut.build_lut(assignments, labels, n_clusters, n_classes) → 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.

Returns ndarray, (K, C), float64 — each row is a valid probability distribution summing to 1.0.

---

reorder_lut

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

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

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

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

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

---

lookup_single_zone

bmc.core.inference.lookup_single_zone(vectors, seeds, lut, metric) → ndarray, (M, C)

Hard LUT lookup for one zone: find nearest seed per vector, return its LUT row.

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.

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

---

accumulate_zones

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

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

Convert accumulated distributions to hard predictions via argmax.

---

forward

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), float64.

---

normalize_distributions

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

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)

class bmc.core.interpolators.Interpolator

Abstract methods

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

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

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

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

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

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

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

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

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

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

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

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), float64 — interpolated distributions. Each row sums to 1.0.

---

forward_interpolated

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), float64.

---

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():

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:

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

---

bmc.backend — Device Backend

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:

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

---

resolve_device

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

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) float32	L2 distance via expand trick in float64.
pairwise_cosine	(A, B, eps=1e-8) → (M, K) float32	Cosine distance.
pairwise_frobenius	(A, B) → (M, K) float32	Frobenius distance (flattens inputs).
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) float32	Runs in float32 on GPU.
pairwise_cosine	(A, B, eps=1e-8) → (M, K) float32	Normalize + matmul on GPU.
pairwise_frobenius	(A, B) → (M, K) float32	Identical to pairwise_l2 (flattened).
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:

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.

---

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.