"""atom_hifi/training/cluster_base.py — backend-agnostic clustering primitives.
The pieces every cluster backend shares, independent of which clusterer runs:
the PCA embedding (:class:`PCAEmbedder`), the deployable classifier
(:class:`NearestNeighborClassifier`), the per-element result contract
(:class:`ClusterFit`), and the trailing-noise-slot promotion machinery
(issue #63). Backend implementations live in
:mod:`atom_hifi.training.cluster_backends`; this module must not import them
(the dependency runs one way, backends -> base).
Split out of ``cluster_backends.py`` in issue #69; the native classifier landed
in issue #76.
"""
import contextlib
import copy as _copy
import logging
import os
import pickle
import tempfile
from abc import ABC, abstractmethod
from dataclasses import dataclass
import numpy as np
from scipy.spatial.distance import cdist
import decaf
from atom_hifi.training.lambda_open import check_lambda
logger = logging.getLogger(__name__)
#: HDBSCAN noise fraction above which the softness/superionic warning fires
#: (cf. LLZO Li noise-frac 0.42 in the §39 three-way study).
NOISE_FRAC_WARN = 0.3
#: distance-matrix element budget per prediction chunk (a test shrinks this)
CHUNK_MAX_EL = int(5e7)
# ── The three abstractions (ADR-0011, issue #69) ──────────────────────────────
[docs]
class Embedder(ABC):
"""Descriptors → the metric space the clusterer and δ both live in.
**Closed on purpose — do not add a second implementation.** δ = c·ε_noise,
the ``|S|`` targets and the QHO σ table are all *distances in decaf's PCA
metric* (broken-stick dim, uncentred projection). A different embedder is a
change of units, not a swap, and no artifact records which embedder produced
it — caches tag the backend only. Embedder-tagging is the prerequisite for
opening this up, not a detail.
"""
@property
@abstractmethod
def t_embedded(self):
"""(N, dim) embedded training rows."""
@property
@abstractmethod
def embed_op(self):
"""(D, dim) projection operator; ``embed`` is a matmul with it."""
[docs]
@abstractmethod
def embed(self, desc):
"""Project raw descriptor rows into the embedding."""
[docs]
class Classifier(ABC):
"""Assigns novel embedded rows to fitted groups, abstaining when far.
Internal abstraction. The *mechanism* is backend-agnostic and, since issue
#76, actually written that way — see :class:`NearestNeighborClassifier`.
Loop and patient OOD consume either supported mechanism. Concretely:
- decaf sets ``ClusterFit.abstain_radius`` -> the 1-NN-over-training-points
rule (:class:`NearestNeighborClassifier`).
- kmeans sets ``ClusterFit.abstain_lambda`` -> the nearest-centroid rule
priced by the DP-means cluster-opening price (gitlab#80,
:class:`NearestCentroidClassifier`). kmeans' granularity is still a
cluster count, but the FITTED PARTITION has a resolution, and
``lambda_open`` prices it — that is what changed.
- hdbscan stores frozen density hierarchy prediction state. It uses
mutual-reachability attachment, not a radius or a centroid price.
``ClusterBackend.has_classifier`` is True for all three built-ins. Kmeans
``lambda_open`` tracks local partition spacing, so when k grows with the
corpus its OOD rate can stay approximately scale-free and the loop may
consume the full library. Gitlab#81 accepts that existing termination
behavior; :meth:`HiFiWorkflow.run_loop` warns once instead of changing the
predictor or inventing an unvalidated stop threshold.
"""
[docs]
@abstractmethod
def classify(self, embedded, **kwargs):
"""Return decaf's ``(gID_distinct, gID_label_all)`` pair.
Both lists have one slot per selection group plus a trailing slot.
``gID_distinct[-1]`` holds rejected rows. DECAF and KMeans provide
unconditional assignments in ``gID_label_all`` (empty trailing slot).
HDBSCAN has no unconditional assignment: rejected rows also occupy
its second list's trailing slot. Promoted noise slots stay empty.
"""
[docs]
class ClusterBackend(ABC):
"""The clusterer, and the one public extension point.
Capabilities are **class attributes**, not name literals: code asks what a
backend can do, never what it is called. ``name`` is the exception — it is
the cache tag, and ``state_io`` enforces equality on it across runs.
"""
#: cache tag — must be stable across runs (state_io enforces equality)
name: str
#: 'distance' (δ is a radius) | 'count' (δ is a cluster count k)
granularity_kind: str
#: fit yields a deployable classifier → OOD / patient / loop
has_classifier: bool = False
#: fit can leave a trailing noise slot needing ε_noise promotion (#63)
emits_noise: bool = False
#: fit yields its own bandwidth estimate (decaf's HaarWavelet δ_auto)
has_auto_bandwidth: bool = False
[docs]
@abstractmethod
def fit(self, desc_el, granularity, *, t_embedded=None, **kw):
"""Cluster one element; return an unpromoted :class:`ClusterFit`.
Takes raw ``desc_el`` rather than a ready embedding because decaf builds
its own embedding internally — ``inject_precomputed`` needs the decaf
object itself. Backends that cluster in the shared embedding should
call :meth:`shared_embedder` and use its ``t_embedded``; that is the
one-liner a custom backend needs.
"""
[docs]
@classmethod
def unsupported_reason(cls, granularity=None):
"""Why this backend cannot run *in this environment*, or None if it can.
The capabilities above are static facts about the algorithm; this one is
not. It reports a dependency fault — a broken upstream version, a
missing optional package — so it has to be asked at runtime rather than
read off the class.
Two questions, one method:
- ``granularity=None`` asks the general one — *can this backend be
relied on here at all?* Scan entry points and the test skip marker
use this, because a scan sweeps δ and will reach the bad values.
- passing a δ asks about **that one fit**, which lets a backend admit
the δ values it can still serve. ``fit_cluster_backend`` uses this,
so a dependency fault that only bites part of the δ range does not
take the whole backend down with it.
"""
return None
[docs]
@staticmethod
def shared_embedder(desc_el, t_embedded=None):
"""The shared PCA embedding: reuse the SHM-served one, else fit it."""
if t_embedded is not None:
return PCAEmbedder(t_embedded, desc_for_lazy=desc_el)
return PCAEmbedder.fit(desc_el)
# ── Persistence (issue #87) ───────────────────────────────────────────────────
[docs]
class PickleIO:
"""Pickle persistence with a **declared** footprint, shared by all backends.
A mixin, not a fourth abstraction: it declares no abstract method, so it
stays a plain base class and leaves its subclasses' metaclass alone.
Backend-neutral by construction: pickle serialises whatever the object
holds, so one implementation here covers decaf, kmeans, hdbscan and any
registered custom backend with no per-backend code — the reason it lives
beside the abstractions rather than on ``Judge``.
What it adds over a bare ``pickle.dump`` is the footprint contract. A fit
artifact keeps references to arrays that are *inputs*, not results (raw
descriptors above all), and pickling them silently turns a 10 MB object into
a 300 MB file. Subclasses list those attributes in ``_PICKLE_DROP``; they
come back as ``None`` and the loader is expected to say so. Saves stamp the
schema and package version so a stale file is diagnosable rather than
mysterious.
The footprint applies to **serialization only**. ``__getstate__`` also
drives ``copy.copy``/``copy.deepcopy``, so without the overrides below a
caller who merely copied a fitted object would silently lose the dropped
attributes on a live object no one asked to persist. Those overrides assume
``__dict__``-based subclasses, which is what ``Judge`` is.
Serialization never *computes*: these objects also travel over joblib IPC,
where a hook that materialised a lazy field would cost real throughput.
"""
#: Bump when a load of an older file would be wrong rather than merely thin.
#: 1 — first stamped format; ``_PICKLE_DROP`` attributes restore as None.
PICKLE_SCHEMA = 1
#: Attribute names nulled on the way out (rebuildable inputs, not results).
_PICKLE_DROP = ()
def __getstate__(self):
from atom_hifi import __version__
state = dict(self.__dict__)
for name in self._PICKLE_DROP:
if name in state:
state[name] = None
state['_pickle_schema'] = self.PICKLE_SCHEMA
state['_atom_hifi_version'] = __version__
return state
def __copy__(self):
clone = self.__class__.__new__(self.__class__)
clone.__dict__.update(self.__dict__)
return clone
def __deepcopy__(self, memo):
clone = self.__class__.__new__(self.__class__)
memo[id(self)] = clone
for name, value in self.__dict__.items():
clone.__dict__[name] = _copy.deepcopy(value, memo)
return clone
[docs]
def save(self, path):
"""Pickle this object to ``path`` (footprint per ``_PICKLE_DROP``).
Written to a temporary sibling and moved into place with
:func:`os.replace`, because ``open(path, 'wb')`` truncates *first*: a
save that failed halfway — the very case the workflow's save site
anticipates for an unpicklable custom backend — would otherwise destroy
the artifact a previous run wrote and leave a truncated file that fails
at load time instead of here.
"""
path = os.fspath(path)
directory = os.path.dirname(path) or '.'
fd, tmp = tempfile.mkstemp(dir=directory,
prefix=f'.{os.path.basename(path)}.',
suffix='.tmp')
try:
with os.fdopen(fd, 'wb') as f:
pickle.dump(self, f, protocol=pickle.HIGHEST_PROTOCOL)
os.replace(tmp, path)
except BaseException:
with contextlib.suppress(OSError):
os.unlink(tmp)
raise
[docs]
@classmethod
def load(cls, path):
"""Unpickle an object saved by :meth:`save`.
A file older than the schema is loaded, not rejected: the whole point of
the artifact is to survive the run that wrote it, and every migration so
far has been additive. A mismatch warns so the reader knows which
version's semantics they are looking at.
"""
with open(path, 'rb') as f:
obj = pickle.load(f)
schema = getattr(obj, '_pickle_schema', None)
if schema is None:
logger.debug(f'{path}: unstamped pickle (written before schema '
f'{cls.PICKLE_SCHEMA}).')
elif schema != cls.PICKLE_SCHEMA:
logger.warning(
f'{path} was written with pickle schema {schema}; this '
f'atom-hifi expects {cls.PICKLE_SCHEMA} '
f'(file written by atom-hifi '
f'{getattr(obj, "_atom_hifi_version", "unknown")}).')
return obj
# ── Implementations of the abstractions ───────────────────────────────────────
[docs]
class PCAEmbedder(Embedder):
"""decaf PCA embedding, computed once per element and shared by all backends.
``embed`` is a pure matmul with ``embed_op`` (decaf's ``embed_PCA``), so a
precomputed (t_embedded, embed_op) pair reproduces decaf's embedding exactly.
On the SHM fast path only ``t_embedded`` is available; ``embed_op`` is then
recomputed lazily from the raw descriptors on first access (scan workers
never touch it, so the extra PCA is never paid in the hot loop).
"""
def __init__(self, t_embedded, embed_op=None, desc_for_lazy=None):
# Backing field, not a plain attribute: Embedder declares t_embedded
# abstract, and an instance attribute set in __init__ does not clear
# __abstractmethods__ — the class would be uninstantiable.
self._t_embedded = np.asarray(t_embedded)
self._embed_op = None if embed_op is None else np.asarray(embed_op)
self._desc_for_lazy = desc_for_lazy
[docs]
@classmethod
def fit(cls, desc_array):
"""Fit decaf's PCA on ``desc_array`` (no clustering)."""
ec = decaf.embed_cluster(np.asarray(desc_array),
embed_str='PCA', cluster_str='MSC')
ec.get_embed(dim=0)
return cls(np.asarray(ec.t_embedded), np.asarray(ec.embed_op))
@property
def t_embedded(self):
return self._t_embedded
@t_embedded.setter
def t_embedded(self, value):
# Plain assignment was valid before #69 (it was a bare attribute).
self._t_embedded = value
# No __getstate__ here on purpose. Dropping `_desc_for_lazy` — the raw
# descriptor reference — looks like the obvious footprint win, but it is
# unreachable: it is only ever set when `_embed_op` is None
# (`shared_embedder`), and the first `embed_op` access clears it, so the two
# are never live together. A hook that dropped it would save nothing, while
# hand-rolling the ``(__dict__, slots)`` protocol it must preserve loses
# name-mangled private slots that the default protocol handles. The real
# footprint is stated once, on `Judge._PICKLE_DROP`. (Materialising
# `embed_op` during serialization to force the drop is not an option: scan
# workers pickle embedders over joblib IPC, and an eager PCA per hand-off is
# a throughput regression.)
def __setstate__(self, state):
# A PCAEmbedder pickled before #69 stored `t_embedded` directly. The
# property reads `_t_embedded`, so without this the unpickled object
# loads fine and then raises AttributeError on first access — the
# failure appears at use, far from the load. Migrate on the way in.
# The default reduce protocol emits ``(__dict__, slot_state)`` for a
# subclass with slots. Preserve that protocol instead of assuming all
# state is one mapping.
if isinstance(state, tuple):
dict_state, slot_state = state
else:
dict_state, slot_state = state, None
dict_state = dict(dict_state or {})
if '_t_embedded' not in dict_state and 't_embedded' in dict_state:
dict_state['_t_embedded'] = dict_state['t_embedded']
# If both spellings exist, the current backing field is authoritative.
dict_state.pop('t_embedded', None)
self.__dict__.update(dict_state)
for name, value in (slot_state or {}).items():
setattr(self, name, value)
@property
def embed_op(self):
if self._embed_op is None:
if self._desc_for_lazy is None:
raise RuntimeError(
'PCAEmbedder: embed_op unavailable (built from a bare '
'precomputed embedding without descriptors).')
self._embed_op = _pca_embed_op(self._desc_for_lazy)
self._desc_for_lazy = None
return self._embed_op
[docs]
def embed(self, desc):
return np.asarray(desc) @ self.embed_op
[docs]
class NearestNeighborClassifier(Classifier):
"""1-NN label transfer over training points + radius abstention.
The native implementation ADR-0011's abstraction table names. It replaces
the delegation to decaf, and is a **verbatim** port of ``MSC_classify`` on
the one path atom_hifi uses (``classify_mode='TrainingData'``,
``label_all='Both'``; ``decaf.py:1510-1554``) — not a reimplementation of
the idea. ADR-0011 left this open precisely because a loose port moves OOD
counts silently, so three conventions are copied rather than improved:
- **``argsort(axis=1)[:, 0]``, not ``argmin``** (``decaf.py:1508/1518``).
numpy's default sort is not stable, so on an exact distance tie the two
*may* resolve to different training atoms — and hence different groups.
Whether they actually differ is numpy-version-dependent (2.2 diverges on
an all-tied row, 2.4/2.5 do not), which is precisely why the call is
copied rather than reasoned about: copying makes this class agree with
decaf on every numpy, including future ones where the tie-break changes
again. Exact ties are not hypothetical in a symmetric crystal.
``argmin`` would also be cheaper in memory; that trade is deferred to
gitlab#77, which measures tie frequency first.
- **``<=``, not ``<``** for the radius test (``decaf.py:1542``).
- The abstention test is over the **training set**, not the centers, and is
therefore independent of how the partition was drawn (``:1540-1542``).
That independence is exactly why this class is backend-agnostic.
Inputs are ``(t_embedded, labels, n_groups, radius)``. Note ``labels`` is
the fit's **promoted** label array, where decaf reads ``t_gID`` through
``reformat_groupID``. The two agree on the production path — ``MeanShift``
runs with ``cluster_all=True`` (``decaf.py:1434``) so no atom is ever ``-1``
and the trailing slot stays empty — but they would not agree after issue-#63
promotion renumbered a backend that does emit noise. Using the promoted
labels is what keeps a promoted fit self-consistent.
"""
def __init__(self, t_embedded, labels, n_groups, radius):
radius = float(radius)
if not np.isfinite(radius) or radius < 0:
# Pre-emptive, NOT a live bug: today the only setter of
# `abstain_radius` is DecafBackend.fit, from a fitted bandwidth, so
# no production path can reach here with a NaN. It becomes live the
# moment gitlab#77 wires a count-kind radius, because
# `HiFiWorkflow._fit_iteration` writes eps_noise = nan for a
# count-kind run that has no ε_noise, and `d <= nan` is False for
# every d — 100 % of atoms marked out-of-distribution, no error
# anywhere. Refused here, where the radius is consumed, rather than
# at each future call site.
raise ValueError(
f'abstention radius must be finite and non-negative; got '
f'{radius!r}. A NaN radius marks every atom out-of-distribution.')
self.t_embedded = np.asarray(t_embedded)
self.labels = np.asarray(labels)
self.n_groups = int(n_groups)
self.radius = radius
[docs]
def classify(self, embedded, dist_in=None, **kwargs):
"""Return decaf's ``(gID_distinct, gID_label_all)`` pair.
``dist_in`` is a precomputed ``(N_new, N_train)`` distance matrix. Both
"absent" spellings are accepted: ``None``, and decaf's own empty-array
sentinel (``decaf.py:1504/1511``), so a call ported straight off a decaf
one does not land on a shape error deep inside ``argsort``. A *present*
matrix of the wrong shape is rejected, which decaf does not do — see the
branch below for why that cannot move a number.
The multiscale arguments are deliberately NOT accepted, because nothing
in atom_hifi passes them and quietly ignoring a multiscale request would
answer a different question than the caller asked.
"""
if kwargs:
raise TypeError(
f'NearestNeighborClassifier.classify got unsupported argument(s) '
f'{sorted(kwargs)}; only dist_in is supported (decaf\'s '
'multibandwidth/multiscales paths are not ported).')
X = np.asarray(embedded)
if dist_in is None or np.asarray(dist_in).shape[0] == 0:
dist = cdist(X, self.t_embedded)
else:
dist = np.asarray(dist_in)
# The one place this class is deliberately STRICTER than decaf.
# decaf validates nothing here and both likely mistakes are silent,
# not loud: too few columns scans only part of the training set for
# abstention — the identical defect to a truncated scan, which is a
# mutation this module's tests exist to catch — and a mismatched row
# count returns answers for rows the caller never asked about.
# Refusing them cannot move a number, because no production caller
# passes dist_in at all (runner._score_remaining_ood,
# hifi_workflow._classify_patient_ood, HiFiWorkflow.run_loop and
# scripts/predict_ood.py::_classify_el all call classify(embedded));
# this argument exists only so a call ported off a decaf one keeps
# working.
n_train = len(self.t_embedded)
if dist.ndim != 2 or dist.shape[1] != n_train:
raise ValueError(
f'dist_in must be a 2-D (N_new, N_train) distance matrix '
f'with N_train={n_train}; got shape {dist.shape}.')
if X.ndim == 2 and dist.shape[0] != X.shape[0]:
raise ValueError(
f'dist_in has {dist.shape[0]} rows but embedded has '
f'{X.shape[0]}; they must describe the same points.')
nearest = dist.argsort(axis=1)[:, 0]
transferred = self.labels[nearest]
groups = [np.where(transferred == g)[0] for g in range(self.n_groups)]
# decaf builds the label-all variant by list concatenation *before*
# rebinding n_gID (decaf.py:1528-1530/:1545), so it keeps the
# PRE-setdiff group arrays: the two returned lists share no array, and
# gID_label_all still contains the rows that abstained. Preserved here
# — which is why `groups` is captured first and setdiff'd into a NEW
# list below rather than filtered in place. The trailing slot is
# dtype=int like decaf's, not the float np.empty((0)) default.
gID_label_all = groups + [np.empty((0), dtype=int)]
# decaf transposes to (N_train, N_new) and tests each column; over the
# untransposed matrix that is exactly "any training point within radius".
distinct = np.where(~(dist <= self.radius).any(axis=1))[0]
gID_distinct = [np.setdiff1d(item, distinct) for item in groups] + [distinct]
return gID_distinct, gID_label_all
def _chunk_slices(n_rows, k):
"""Row ranges of a (n_rows, k) distance matrix within the element budget.
The budget is per DISTANCE MATRIX, not peak memory, and it has one honest
hole: ``max(1, CHUNK_MAX_EL // k)`` cannot go below one row, so when a
SINGLE row already exceeds the budget (k > CHUNK_MAX_EL) the yielded chunk
is k elements wide regardless. That needs more than 5e7 centroids for one
element, so it is unreachable in practice and is documented rather than
guarded — a guard would have to either raise on a legal input or silently
split along the centroid axis, which would change the argmin. Callers
asking for confidence also pay a second chunk-sized temporary inside
``np.partition``.
"""
if n_rows * k <= CHUNK_MAX_EL:
yield slice(0, n_rows)
return
chunk = max(1, CHUNK_MAX_EL // k)
for start in range(0, n_rows, chunk):
yield slice(start, min(start + chunk, n_rows))
[docs]
class NearestCentroidClassifier(Classifier):
"""Nearest-centroid label transfer + DP-means opening-price abstention.
The kmeans counterpart of :class:`NearestNeighborClassifier`. The
assignment rule IS kmeans' own definition of a label — ``argmin_j
||x - mu_j||`` — which is the point: the 1-NN-over-training-points rule
assigns a DIFFERENT cluster on 10-41% of the queries both rules resolve
(gitlab#80 Stage A, measured on 7 real element arms across 3 systems),
so a kmeans fit cannot be answered through NearestNeighborClassifier.
Abstention uses the DP-means cluster-opening price: q(x) = d1**2 /
lambda_open, with q == 1 EXACTLY in distribution (STRICTLY greater
means OOD). Note the units: `lambda_open` is in SQUARED-distance
units — d1**2, not d1, is what is compared to it.
"""
def __init__(self, centers, n_groups, lambda_open):
# check_lambda FIRST, before any point can be classified: a
# NaN/inf/negative/zero lambda fails hard here, and each failure's
# message states the consequence of classifying under it.
self.lambda_open = check_lambda(
lambda_open, context='NearestCentroidClassifier')
centers = np.asarray(centers, dtype=float)
if centers.ndim != 2:
raise ValueError(
f'centers must be a 2-D (n_groups, emb_dim) array; got shape '
f'{centers.shape}.')
if len(centers) != int(n_groups):
raise ValueError(
f'got {len(centers)} centers for n_groups={int(n_groups)}; '
'group ids index centers directly, so a mismatch would '
'silently mislabel every row.')
if not np.isfinite(centers).all():
bad = np.flatnonzero(~np.isfinite(centers).all(axis=1))
raise ValueError(
f'centers rows {bad.tolist()[:8]} are not finite. A NaN in a '
'centre makes every cdist distance to it NaN, argmin then '
'selects that column (NaN comparisons are all False, so the '
'first column wins), and the abstention test q > 1 is False '
'for q = NaN — so EVERY query would be labelled with the '
'broken group and reported as in-distribution. Non-finite '
'geometry must be traced, never classified.')
self.centers = centers
self.n_groups = int(n_groups)
[docs]
def predict(self, embedded, return_confidence=False):
"""Centroid rule: assign nearest centroid, abstain where q > 1.
q(x) = d1**2 / lambda_open, with the boundary q == 1 exactly IN
distribution (STRICTLY greater means OOD). d2 is the second-smallest
distance (inf and confidence 1.0 when k == 1); confidence is
``1 - d1/d2``, 0 where d1 == d2 == 0. Distances are computed in row
chunks (each chunk's (rows, k) matrix under CHUNK_MAX_EL elements) so
memory stays bounded; d1 and d2 come from the SAME distance pass.
Returns a dict with ``label`` (-1 where OOD), ``assigned`` (the raw
nearest-centroid index, never -1), ``d1``, ``d2``, ``q``,
``confidence``, ``ood`` and ``ood_fraction``. ``d2`` and
``confidence`` are None unless ``return_confidence``.
"""
X = np.asarray(embedded)
if len(X) and not np.isfinite(X).all():
bad = np.flatnonzero(~np.isfinite(np.atleast_2d(X)).all(axis=1))
raise ValueError(
f'{len(bad)} query rows are not finite (first: '
f'{bad.tolist()[:8]}). d1 would be NaN for them, and the OOD '
'test q > 1 is False for NaN, so an invalid row would be '
'reported as confidently represented by group 0 instead of '
'being rejected — the exact opposite of the truth. Trace the '
'embedding before classifying.')
k = self.n_groups
assigned = np.empty(len(X), dtype=np.int64)
# SQUARED distances throughout the pass. `lambda_open` is in squared
# units, so q must be formed from a squared distance that was never
# square-rooted: `cdist(...)` takes a sqrt and `d1 * d1` puts it back,
# and that round-trip is not the identity. Measured on a real
# calibration state (an occupancy-2 cluster whose charged centre sits
# on one member, where lambda == r100**2 EXACTLY because the m/(m-1)
# slack is consumed by the drift term): the direct squared residual is
# 11.436327734832187 and the round-trip gives ...188, one ULP higher,
# so q came out 1.0000000000000002 and the classifier ABSTAINED ON ONE
# OF ITS OWN TRAINING ATOMS — the single contract gitlab#80 exists to
# provide. 'sqeuclidean' returns the sum of squared differences
# directly and reproduces the calibration's own arithmetic on that
# state. Ordering is unchanged (t -> t**2 is monotone on t >= 0), so
# argmin and the second-smallest are the same rows as before.
d1_sq = np.empty(len(X))
d2_sq = np.empty(len(X)) if return_confidence else None
for sl in _chunk_slices(len(X), k):
D = cdist(X[sl], self.centers, 'sqeuclidean')
# argmin, NOT argsort(axis=1)[:, 0]. The argsort convention in
# NearestNeighborClassifier exists solely to reproduce the
# third-party mean-shift implementation bit-for-bit on exact
# distance ties (see its constructor docstring); there is no decaf
# counterpart for the centroid rule, and np.argmin
# deterministically returns the FIRST minimum on a tie, which is
# the stronger guarantee.
idx = D.argmin(axis=1)
assigned[sl] = idx
d1_sq[sl] = D[np.arange(sl.stop - sl.start), idx]
if return_confidence:
# np.partition(D, 1) puts the second-smallest at column 1.
d2_sq[sl] = (np.partition(D, 1, axis=1)[:, 1] if k >= 2
else np.full(len(idx), np.inf))
# The public d1/d2 stay DISTANCES — that is the documented contract and
# what `confidence = 1 - d1/d2` is defined on. Only q is formed from
# the squared values, and the sqrt here is over n values rather than
# the n*k matrix cdist used to sqrt, so this is also cheaper.
d1 = np.sqrt(d1_sq)
d2 = np.sqrt(d2_sq) if return_confidence else None
q = d1_sq / self.lambda_open
ood = q > 1.0 # strictly greater; q == 1 exactly is in distribution
label = assigned.copy()
label[ood] = -1
confidence = None
if return_confidence:
# d2 == 0 implies d1 == 0 (d2 >= d1 always) — the 0/0 case is a
# query exactly on a centre with a coincident duplicate centre.
# np.divide with `where=` rather than np.where(cond, 0.0, d1/d2):
# the latter evaluates the division for EVERY row before selecting,
# so the 0/0 rows still execute and emit
# "RuntimeWarning: invalid value encountered in divide" — a warning
# about a case this line already handles correctly.
confidence = np.zeros_like(d1)
np.divide(d1, d2, out=confidence, where=d2 != 0)
confidence = np.where(d2 == 0, 0.0, 1.0 - confidence)
return dict(
label=label,
# The raw argmin, WITHOUT the -1 abstention marker written over it.
# `classify` needs it: decaf's gID_label_all labels every row
# including the abstained ones, so overwriting in place would force
# a second distance pass over the OOD rows to recover what this one
# already computed. Kept as its own key rather than by making
# `label` mean two things.
assigned=assigned,
d1=d1,
d2=d2,
q=q,
confidence=confidence,
ood=ood,
# 0.0 for an empty query set — mean over an empty array is NaN.
ood_fraction=float(ood.mean()) if len(ood) else 0.0,
)
[docs]
def classify(self, embedded, **kwargs):
"""Return decaf's ``(gID_distinct, gID_label_all)`` pair.
``dist_in`` is deliberately NOT accepted: it is a ``(N_new, N_train)``
distance matrix over TRAINING POINTS, which is meaningless for a
centroid rule — there is no training-point distance anywhere in this
classifier.
"""
if kwargs:
raise TypeError(
f'NearestCentroidClassifier.classify got unsupported '
f'argument(s) {sorted(kwargs)}; in particular dist_in is NOT '
'accepted: it is a (N_new, N_train) distance matrix over '
'TRAINING POINTS, which is meaningless for a centroid rule.')
X = np.asarray(embedded)
# Confidence costs an extra O(M k) partition pass and must not be
# paid by the abstention path.
result = self.predict(X, return_confidence=False)
ood = result['ood']
# `assigned`, not `label`: decaf's gID_label_all labels EVERY row by
# its group, abstained rows included (decaf.py:1528-1530/:1545), and
# `label` carries the -1 marker over exactly those rows.
assigned = result['assigned']
groups = [np.where(assigned == g)[0] for g in range(self.n_groups)]
# Same construction as NearestNeighborClassifier.classify: the
# PRE-setdiff group arrays go into gID_label_all, and gID_distinct
# setdiffs into a NEW list, so the two returned lists share no array.
gID_label_all = groups + [np.empty((0), dtype=int)]
distinct = np.flatnonzero(ood)
gID_distinct = [np.setdiff1d(g, distinct) for g in groups] + [distinct]
return gID_distinct, gID_label_all
[docs]
@dataclass
class ClusterFit:
"""Per-element clustering result — the backend-agnostic contract.
Fits returned by :func:`fit_cluster_backend` are self-consistent: every
label is in ``[0, n_groups)``, ``centers`` and ``noise_derived`` have
``n_groups`` rows, and the trailing noise/distinct slot is empty.
"""
labels: np.ndarray # (N_el,) int
n_groups: int # density-supported + promoted noise groups
centers: np.ndarray # (n_groups, emb_dim)
granularity: float # δ or k actually used
backend: str
embedder: Embedder
#: Optional :class:`Classifier` **override**. Since issue #76 no backend
#: sets it: :meth:`classify` builds a :class:`NearestNeighborClassifier`
#: from ``abstain_radius`` instead, so nothing stateful rides on the fit.
#: It stays honoured first, which keeps an explicitly-supplied
#: implementation — and any pickle written before #76 — working unchanged.
#:
#: Renamed from ``decaf_cluster`` in issue #69 — deliberately a NEW name,
#: not a repurposed one. A pickle written between #68 and #69 still
#: unpickles (a dataclass restores __dict__ wholesale, so the stale
#: ``decaf_cluster`` key simply rides along); had this field kept that name,
#: the stale value would satisfy the `is not None` test and the fit would
#: classify through a stale decaf object with no error at all. Under the
#: new name such a fit raises, which is the honest outcome.
classifier: object = None
# (n_groups,) bool — True for a group created by promoting a deduped noise
# representative (issue #63). All-False when no promotion occurred.
noise_derived: np.ndarray = None
# Noise fraction measured BEFORE promotion. Afterwards no label is −1, so
# the raw computation would read 0 and silently disable the softness flag
# (cf. LLZO 0.42).
pre_promotion_noise_frac: float = None
#: Radius for the 1-NN abstention test, or None when this backend has no
#: deployable classifier. It carries the *policy*, which is why it is a
#: field and not derived from ``granularity``: for decaf the two coincide
#: (δ IS decaf's radius), but kmeans' granularity is a cluster count with no
#: distance meaning. ``DecafBackend`` sets this, kmeans sets
#: ``abstain_lambda`` (gitlab#80), and hdbscan stores density state instead.
abstain_radius: float = None
#: DP-means cluster-opening price in SQUARED-distance units (gitlab#80),
#: set by ``KMeansBackend.fit`` only. It is a SEPARATE field from
#: ``abstain_radius`` and not a unit conversion of it, because the two
#: select *different classifiers* — nearest-centroid vs 1-NN-over-training
#: points — and Stage A measured those two disagreeing on 10-41% of
#: jointly-resolved queries. ``ClusterFit.classify`` raises if both are
#: set, rather than silently preferring one.
abstain_lambda: float = None
#: Provenance dict returned by ``lambda_local`` (gitlab#80), persisted so
#: a fit's abstention threshold can be audited after the fact.
abstain_calibration: dict = None
#: Versioned frozen density state, built from sklearn's fitted hierarchy.
#: Separate from promotion: its labels remain backend-local density ids.
hdbscan_prediction: object = None
@property
def n_density_groups(self):
"""Coverable groups produced by clustering, excluding noise promotion."""
promoted = (int(np.count_nonzero(self.noise_derived))
if self.noise_derived is not None else 0)
return self.n_groups - promoted
@property
def n_deployable_groups(self):
"""Number of targets supported by this saved fit's prediction state."""
if self.classifier is not None:
return self.n_groups
if self.backend == 'hdbscan':
state = self.hdbscan_prediction
return 0 if state is None else state.n_density_groups
if self.abstain_radius is not None or self.abstain_lambda is not None:
return self.n_groups
return 0
# The decaf-compat accessors are gone (issue #69): groupnum, cluster_centers,
# bandwidth, t_embedded, embed_op, embed and t_gID. Use the native fields —
# n_groups, centers, granularity, embedder.* — and index `labels` directly
# instead of t_gID, which rebuilt every group slice on every access.
# `bandwidth` was a lie on the count-kind path: it returned a cluster count.
[docs]
def classify(self, embedded, **kwargs):
"""Assign novel embedded rows to this fit's groups, abstaining when far.
The classifier is built here rather than stored on the fit — whether
:class:`NearestNeighborClassifier` (from ``abstain_radius``) or
:class:`NearestCentroidClassifier` (from ``abstain_lambda``). It holds
only references plus a float, so construction is O(1) and there is
nothing to amortise — whereas a cached instance would put a second
reference to ``t_embedded`` on an object ``Judge.save`` pickles. This
is also what closes ADR-0011's open pickle-footprint item: no live decaf
object needs retaining on a fit any more.
``classifier`` is still honoured first, so an explicitly-supplied
implementation (and any pickle written before this change) keeps working.
"""
if self.classifier is not None:
return self.classifier.classify(embedded, **kwargs)
if self.backend == 'hdbscan':
if self.hdbscan_prediction is None:
raise RuntimeError(
'This legacy hdbscan fit has no frozen prediction state '
'and no deployable classifier. Explicitly fit a new Judge '
'from the original reference descriptors and save it with '
'this version; loading never refits (issue #97).')
return self.hdbscan_prediction.classify(
embedded, n_groups=self.n_groups, **kwargs)
if self.abstain_lambda is not None and self.abstain_radius is not None:
# They encode two different classifiers — nearest-centroid vs
# 1-NN-over-training-points, which disagree on 10-41% of
# jointly-resolved queries (gitlab#80) — so silently preferring
# one would answer a different question than the caller asked.
raise ValueError(
f'{self.backend!r} fit sets BOTH abstain_lambda and '
'abstain_radius; they encode two different classifiers '
'(nearest-centroid vs 1-NN-over-training-points, which '
'disagree on 10-41% of jointly-resolved queries, gitlab#80), '
'so silently preferring one would answer a different '
'question than the caller asked.')
if self.abstain_lambda is not None:
return NearestCentroidClassifier(
self.centers, self.n_groups, self.abstain_lambda,
).classify(embedded, **kwargs)
if self.abstain_radius is None:
raise RuntimeError(
f'{self.backend!r} backend has no deployable classifier — it '
'sets neither abstain_radius (1-NN rule, decaf) nor '
'abstain_lambda (nearest-centroid rule, kmeans/gitlab#80), so '
'there is no abstention threshold to classify against '
'(ADR-0009).')
return NearestNeighborClassifier(
self.embedder.t_embedded, self.labels, self.n_groups,
self.abstain_radius,
).classify(embedded, **kwargs)
@property
def noise_frac(self):
"""Fraction of atoms the backend put in the trailing slot.
Reads the pre-promotion value when promotion has run, so the softness
signal survives noise reps becoming real groups (issue #63).
"""
if self.pre_promotion_noise_frac is not None:
return float(self.pre_promotion_noise_frac)
return float((self.labels == -1).sum()) / max(len(self.labels), 1)
[docs]
def dedup_noise_greedy(coords, delta):
"""Greedy single-pass dedup of noise atoms at radius δ (= ε_noise).
Walk the noise atoms; keep one as a DISTINCT rare-environment rep if it is
farther than δ from every already-kept rep, else it is a thermal duplicate
absorbed into that rep. Deterministic (natural order). O(n_noise · n_kept).
Returns
-------
kept_idx : list[int]
Indices (into ``coords``) of the kept distinct reps.
assign : ndarray, shape (len(coords),)
``assign[i]`` is the position **within ``kept_idx``** of the rep that
``coords[i]`` belongs to — identity for kept reps themselves. Used to
absorb thermal duplicates into their rep's group rather than discarding
them (issue #63).
"""
kept_idx, kept = [], None
assign = np.zeros(len(coords), dtype=int)
for i in range(len(coords)):
x = coords[i]
if kept is None:
assign[i] = 0
kept_idx.append(i)
kept = x[None, :]
continue
dists = np.linalg.norm(kept - x, axis=1)
j = int(np.argmin(dists))
if dists[j] > delta:
assign[i] = len(kept_idx)
kept_idx.append(i)
kept = np.vstack([kept, x])
else:
assign[i] = j
return kept_idx, assign
[docs]
def medoids_from_labels(t_embedded, labels, n_groups):
"""Per-cluster medoid (point nearest the cluster mean) in the embedding.
HDBSCAN sets no cluster_centers; medoids give a center where the
surroundings extraction needs one. Noise atoms (label < 0) are ignored.
"""
X = np.asarray(t_embedded)
cents = np.zeros((n_groups, X.shape[1]))
for g in range(n_groups):
m = labels == g
if not m.any():
continue
sub = X[m]
mu = sub.mean(0)
cents[g] = sub[np.argmin(np.linalg.norm(sub - mu, axis=1))]
return cents
def _pca_embed_op(desc_el):
"""Recompute decaf's PCA embed_op from raw descriptors (SHM path helper)."""
ec = decaf.embed_cluster(np.asarray(desc_el), embed_str='PCA',
cluster_str='MSC')
ec.get_embed(dim=0)
return np.asarray(ec.embed_op)