diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 56e0e7f..e3524b4 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -60,6 +60,62 @@ jobs: - name: pytest run: uv run pytest --cov --cov-report=term-missing --cov-report=xml + fast-extra: + name: the accelerated path + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1 + - uses: astral-sh/setup-uv@c771a70e6277c0a99b617c7a806ffedaca235ff9 # v9.0.0 + with: + enable-cache: true + python-version: "3.12" + + # numba is NOT a declared extra of this package, deliberately: uv resolves every extra into + # one lockfile, so declaring it would cap NumPy below 2.5 for the whole project and the test + # matrix would then run against an older NumPy than the package releases against. It is + # installed ad hoc here instead, which is also how a user installs it. + - run: uv sync --extra dev + + # These must PASS, not skip. tests/test_numba_kernel.py is guarded by numba being importable, + # so without this job it would skip everywhere and the second implementation of the hottest + # code in the package would go unchecked. + - name: the differential tests run rather than skip + run: | + uv run --with "numba<0.67" pytest tests/test_numba_kernel.py -v --no-header | tee result.txt + grep -q PASSED result.txt + if grep -q SKIPPED result.txt; then echo "numba tests skipped; numba did not install"; exit 1; fi + + - name: the whole suite still passes with numba present + run: uv run --with "numba<0.67" pytest -q + + no-extras: + name: numpy only, no extras + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1 + - uses: astral-sh/setup-uv@c771a70e6277c0a99b617c7a806ffedaca235ff9 # v9.0.0 + with: + enable-cache: true + python-version: "3.12" + + # The claim the package leads with, checked rather than asserted: a plain install pulls NumPy + # and nothing else, and it resolves the *current* NumPy rather than the one numba caps us to. + - name: a plain install is numpy and nothing else + run: | + uv venv /tmp/plain + VIRTUAL_ENV=/tmp/plain uv pip install . + VIRTUAL_ENV=/tmp/plain uv pip list --format=freeze | grep -v '^python-som==' > installed.txt + cat installed.txt + test "$(wc -l < installed.txt)" -eq 1 + grep -q '^numpy==2\.' installed.txt + /tmp/plain/bin/python -c " + import sys, numpy as np, python_som + som = python_som.SOM(x=8, y=8, input_len=3, random_seed=0) + som.train(np.random.default_rng(0).normal(size=(30, 3)), n_iteration=5, mode='batch') + assert 'numba' not in sys.modules + print('numpy-only install OK, numpy', np.__version__) + " + benchmarks: name: benchmarks still run runs-on: ubuntu-latest @@ -96,7 +152,7 @@ jobs: PYTHONPATH: benchmarks run: | uv run python -c " - import bench_update, bench_batch, bench_vs_minisom, bench_vs_sompy + import bench_update, bench_vs_minisom, bench_vs_sompy print('benchmark scripts import OK') " - name: bench_vs_sompy degrades cleanly without its environment diff --git a/CHANGELOG.md b/CHANGELOG.md index a2611dd..8f4fdc8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,43 +5,80 @@ All notable changes to this project are documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). -## [Unreleased] +## [0.7.0] - 2026-07-31 -Nothing in the shipped package changed. This is benchmarking and the tests behind it. +Batch training is 20x to 40x faster and now beats both comparable libraries by a wide margin. The +arithmetic is Kohonen's, reorganised; results move by about 1e-15, which makes this a minor release +rather than a patch. + +### Changed + +- **Batch training is 20x to 40x faster.** Measured against 0.6.1, and against MiniSom and SOMPY + under the fairness protocol from the comparison page: + + | map | 0.6.1 | 0.7.0 | vs MiniSom | vs SOMPY | + | --- | --- | --- | --- | --- | + | 20x20 | 234.2 ms | 12.1 ms | 23.1x faster | 26.2x faster | + | 40x40 | 992.4 ms | 22.7 ms | 30.9x faster | 70.3x faster | + | 60x60 | 2780.2 ms | 54.0 ms | 30.2x faster | 94.4x faster | + | 100x100 | 13972.8 ms | 340.1 ms | | | + + Through 0.6.1 batch training was 1.3x to 1.6x *slower* than MiniSom. Stepwise training is + unchanged and remains about 10% faster than MiniSom's. + + Two changes. Eq. (8) sums over every pair of nodes, and because the neighborhood depends only on + the offset between two nodes that sum is a convolution; both neighborhoods batch training admits + are separable, so it contracts to two matrix products instead of one pass per node. And the + best-matching-unit search expands the Euclidean norm into a single matrix product rather than one + full-grid norm per sample. Kohonen derives Eq. (8) from Eq. (7) on the same grounds, that "the + same addends occur a great number of times" (Section 4.4). + + Peak memory did not rise to pay for it: 3.7 MB against 4.6 MB at 100x100. + +- **Results are not bit-identical to 0.6.1.** The contraction sums the same terms in a different + order, so trained weights differ by 5.9e-16 to 8.5e-16 relative. Below anything a result depends + on, and enough to break an exact-equality check. Pin `python-som==0.6.1` to reproduce an older + figure exactly. ### Added -- **A published comparison against MiniSom and SOMPY**, at - [Comparison with MiniSom and SOMPY](https://andremsouza.github.io/python-som/explanation/comparison-with-som-libraries/). - Both peers implement Kohonen's Eq. (8) and cite the same sources this package works from, so the - comparison is between three implementations of one published algorithm rather than between three - different algorithms. Getting there needed eight controls, because the packages look far more - interchangeable than they are: MiniSom's `train_batch` is stepwise training rather than batch, - only the gaussian is the same function in all three, and SOMPY's - `calculate_quantization_error` returns an elementwise mean absolute error rather than the usual - mean Euclidean distance. - - Measured on batch training, this package is 1.3x to 2.0x faster than SOMPY and 1.3x to 1.6x - **slower** than MiniSom, the latter growing with map size. `batch_update` walks every node in - Python, 108,000 `einsum` calls on a 60x60 map over 30 iterations, where MiniSom's node-side update - is a single vectorised divide. Stepwise training is within 10% of MiniSom's. - -- **`tests/test_minisom_agreement.py`**, which checks the neighborhood functions, the decay, the - best-matching-unit search and both training loops against MiniSom on every run. Trained models - agree to 2.8e-16 relative for stepwise and 1.3e-15 for batch, which is what makes timing the two - against each other meaningful. It also pins the two neighborhoods that deliberately *disagree*, so - neither is later "fixed" into the other. Same reasoning as - `tests/test_linalg_matches_sklearn.py`, which once found a real 5.8% defect in this package. - -- **An asv suite in `asv_benchmarks/`** tracking this package against its own git history, which is - what numpy, scipy, pandas and scikit-learn all use asv for. Sixteen benchmarks over the training - loops, the neighborhood kernel, the matching functions, the SVD and the artifact round trip. No - timing gates anything: CI executes each benchmark once and discards the numbers, because a shared - runner cannot measure anything. - -- **A `bench` extra** holding `asv` and `minisom`, kept out of `dev` because no gating job needs - them. SOMPY is deliberately absent and cannot be added: it uses `np.Inf`, removed in NumPy 2.0, at - class-definition time, so it gets an interpreter of its own that is built by hand. +- **Optional numba acceleration** for the winner search, used automatically when numba is present: + + ```bash + pip install numba + ``` + + Worth 1.0x to 2.4x on top of the above, with bit-identical results, and uneven: where the + neighborhood update dominates it changes little. Deliberately neither a dependency nor an extra. + numba 0.66 requires `numpy<2.5` while this package releases against 2.5, so requiring it would cap + every user's NumPy, and declaring it as an extra caps the lockfile, since uv resolves every extra + together. A plain `pip install python-som` is still NumPy and nothing else, which a CI job checks. + numba is imported on first use, so installing it does not slow `import python_som`. + +- **Two explanation pages**: how batch training is computed, and why linear initialization is an + SVD. Both hold derivations that used to sit in docstrings. + +### Removed + +- **The neighborhood kernel machinery in `python_som._core`**: `gaussian_kernel`, `bubble_kernel`, + `mexican_hat_kernel`, `kernel_view`, `NEIGHBORHOOD_KERNELS`, `resolve_kernel` and `offset_span`. + The axis-matrix contraction replaced them and they had no remaining caller. Private names, so no + supported API changes. + +### Deprecated + +- **`KernelFunction`**, which is in `python_som.__all__` and no longer describes anything the + package produces. It still works and is removed at 1.0.0. + +### Fixed + +- **The winner search is exact for data far from the origin.** Expanding the Euclidean norm is + faster and cancels catastrophically when the models sit far from zero: at an offset of 1e9, + `||w||^2` is around 1e18 while the differences between models are of order 1, and a naive + expansion gives a different node for 499 of 500 samples. Centring both sides on the models' mean + is exact, costs 1%, and removes it at every offset tested up to 1e12. This shipped correct; the + entry is here because it is the same failure mode 0.4.0 fixed in linear initialization, and + because a regression test now pins it. ## [0.6.1] - 2026-07-30 diff --git a/README.md b/README.md index 56b9dc4..739d207 100644 --- a/README.md +++ b/README.md @@ -56,6 +56,8 @@ A full worked example with plots is in [examples/iris.py](https://github.com/and ## Features * NumPy is the only runtime dependency; a fresh install is 69 MB across one package +* Batch training 23x to 31x faster than MiniSom and 26x to 94x faster than SOMPY, measured +* Optional numba acceleration: `pip install numba` and it is used automatically * Stepwise and batch training * Random, random-sampling and linear (PCA) weight initialization * Automatic selection of the map size ratio, from PCA @@ -114,6 +116,12 @@ transitively through this package, depend on them directly, or install `python-s strings are permanent and 1.0.0 will not remove them. If you saw that warning, you can stop migrating. +**0.7.0** makes batch training 20x to 40x faster by reorganising the same arithmetic. Because the +sums happen in a different order, trained weights differ from 0.6.1 by about 1e-15 relative. That is +far below anything a result depends on, and it does break an exact-equality check against a stored +map: pin `python-som==0.6.1` if you need one to match bit for bit. It also removes the neighborhood +kernel helpers from the private `python_som._core`, which had no callers outside the package. + Each change and the passage of Kohonen (2013) behind it is in the [changelog](https://github.com/andremsouza/python-som/blob/master/CHANGELOG.md). @@ -131,6 +139,10 @@ uv run mkdocs serve # docs, locally pre-commit install # optional, run the gates on commit ``` +numba is not a declared extra: uv resolves every extra into one lockfile, so declaring it would cap +NumPy below 2.5 for the whole project. Install it alongside when working on the accelerated path, +with `uv run --with numba pytest tests/test_numba_kernel.py`. + If you use the SonarQube for IDE (SonarLint) VS Code extension, it will also apply Sonar's Python rules locally; the ruff configuration is set up to cover most of the same ground. @@ -140,7 +152,6 @@ Hand-run, never part of the test suite: a timing assertion on shared hardware me ```bash uv run python benchmarks/bench_vs_minisom.py # vs MiniSom, agreement verified before timing -uv run python benchmarks/bench_batch.py # the neighborhood kernel against evaluating per node cd asv_benchmarks && uv run --extra bench asv continuous master HEAD # this package across commits ``` diff --git a/asv_benchmarks/benchmarks/neighborhood.py b/asv_benchmarks/benchmarks/neighborhood.py index 42b9e4b..cbd2295 100644 --- a/asv_benchmarks/benchmarks/neighborhood.py +++ b/asv_benchmarks/benchmarks/neighborhood.py @@ -1,20 +1,22 @@ -"""The neighborhood kernel, which is the optimization 0.4.0 rests on. +"""The axis-matrix contraction, which is how batch training evaluates Eq. (8). -Batch training needs ``h_ji`` for every pair of nodes. Because a neighborhood depends only on the -offset between two nodes, one kernel over every offset serves the whole grid and each node's -neighborhood is a slice of it. Evaluating per node instead was 42% of batch training on a 40x40 map. +Eq. (8) needs ``h_ji`` for every pair of nodes. Because a neighborhood depends only on the offset +between two nodes the sum is a convolution, and both neighborhoods batch training admits are +separable, so it contracts to two matrix products against ``(X, X)`` and ``(Y, Y)`` matrices with no +loop over nodes. -Both paths are still benchmarked, and the per-node one is the point: it is what the kernel -replaced, so keeping it measured is what makes the claim checkable rather than historical. +The per-node path is benchmarked alongside it. That is what the contraction replaced, and keeping it +measured is what makes the claim checkable rather than historical. """ from __future__ import annotations from typing import TYPE_CHECKING -from python_som._core._neighborhood import kernel_view, resolve, resolve_kernel +from python_som._core._neighborhood import axis_matrix, resolve, resolve_axis_profile +from python_som._core._update import batch_update -from .common import RADIUS, SHAPES +from .common import FEATURES, RADIUS, SEED, SHAPES if TYPE_CHECKING: # pragma: no cover from collections.abc import Callable @@ -22,92 +24,111 @@ import numpy as np import numpy.typing as npt -#: Both neighborhoods batch training accepts, plus the signed one only stepwise can use. -NEIGHBORHOODS = ["gaussian", "bubble", "mexican_hat"] +#: Neighborhoods batch training admits, and so the ones with an axis profile. +SEPARABLE = ["gaussian", "bubble"] #: No wrapping, wrapping on one axis, wrapping on both. The cyclic fold is the part of the offset #: machinery most likely to be got wrong, and the part whose cost is least obvious. CYCLIC = [(False, False), (True, False), (True, True)] -class Kernel: - """Building one kernel per iteration.""" +class AxisMatrix: + """Building the two per-axis matrices, once per iteration.""" - params = (SHAPES, NEIGHBORHOODS, CYCLIC) + params = (SHAPES, SEPARABLE, CYCLIC) param_names = ("shape", "neighborhood", "cyclic") - build: Callable[..., npt.NDArray[np.floating]] + profile: Callable[..., npt.NDArray[np.floating]] def setup(self, shape: tuple[int, int], neighborhood: str, cyclic: tuple[bool, bool]) -> None: - """Resolve the kernel builder outside the timed region. + """Resolve the axis profile outside the timed region. :param shape: Grid shape. :param neighborhood: Neighborhood function name. :param cyclic: Whether each axis wraps. """ del shape, cyclic - self.build = resolve_kernel(neighborhood) + self.profile = resolve_axis_profile(neighborhood) def time_build( self, shape: tuple[int, int], neighborhood: str, cyclic: tuple[bool, bool] ) -> None: - """Time building the kernel once. + """Time building both matrices. :param shape: Grid shape. :param neighborhood: Unused. :param cyclic: Whether each axis wraps. """ del neighborhood - self.build(shape, RADIUS, cyclic) + axis_matrix(shape[0], RADIUS, cyclic=cyclic[0], profile=self.profile) + axis_matrix(shape[1], RADIUS, cyclic=cyclic[1], profile=self.profile) def peakmem_build( self, shape: tuple[int, int], neighborhood: str, cyclic: tuple[bool, bool] ) -> None: - """Track the kernel's size, which is the justification for the whole approach. + """Track their size, which is the justification for the approach. - A ``(2X-1, 2Y-1)`` kernel is 198 KB at 80x80 against the 800 MB a full ``(x, y, x, y)`` - tensor would need. A regression here would otherwise be silent. + ``X^2 + Y^2`` floats, against the ``(x, y, x, y)`` tensor the naive form would need, which + reaches 800 MB on a 100x100 map. :param shape: Grid shape. :param neighborhood: Unused. :param cyclic: Whether each axis wraps. """ del neighborhood - self.build(shape, RADIUS, cyclic) + axis_matrix(shape[0], RADIUS, cyclic=cyclic[0], profile=self.profile) + axis_matrix(shape[1], RADIUS, cyclic=cyclic[1], profile=self.profile) -class Slice: - """Taking one node's neighborhood out of a built kernel. +class Contraction: + """One whole Eq. (8) update: both contractions and the guarded divide.""" - Must stay a view rather than a copy. Copying ``(X, Y)`` floats per node would give back most of - what the kernel wins, and this is where that would show up as a trend. - """ + params = (SHAPES, FEATURES) + param_names = ("shape", "n_features") - params = (SHAPES,) - param_names = ("shape",) - - kernel: npt.NDArray[np.floating] - nodes: list[tuple[int, int]] + weights: npt.NDArray[np.floating] + sums: npt.NDArray[np.floating] + counts: npt.NDArray[np.floating] + hx: npt.NDArray[np.floating] + hy: npt.NDArray[np.floating] - def setup(self, shape: tuple[int, int]) -> None: - """Build the kernel and the node list outside the timed region. + def setup(self, shape: tuple[int, int], n_features: int) -> None: + """Build the models, accumulators and axis matrices outside the timed region. :param shape: Grid shape. + :param n_features: Number of features. """ - self.kernel = resolve_kernel("gaussian")(shape, RADIUS, (False, False)) - self.nodes = [(x, y) for x in range(shape[0]) for y in range(shape[1])] + import numpy as np # noqa: PLC0415 asv collects this module without running setup - def time_slice_every_node(self, shape: tuple[int, int]) -> None: - """Time slicing the kernel once per node, which is one batch iteration's worth. + rng = np.random.default_rng(SEED) + self.weights = rng.normal(size=(*shape, n_features)) + self.sums = rng.normal(size=(*shape, n_features)) + self.counts = rng.integers(0, 3, size=shape).astype(float) + profile = resolve_axis_profile("gaussian") + self.hx = axis_matrix(shape[0], RADIUS, cyclic=False, profile=profile) + self.hy = axis_matrix(shape[1], RADIUS, cyclic=False, profile=profile) - :param shape: Grid shape. + def time_update(self, shape: tuple[int, int], n_features: int) -> None: + """Time the contraction. + + :param shape: Unused. + :param n_features: Unused. """ - for node in self.nodes: - kernel_view(self.kernel, shape, node) + del shape, n_features + batch_update(self.weights, self.sums, self.counts, self.hx, self.hy) + + def peakmem_update(self, shape: tuple[int, int], n_features: int) -> None: + """Track what one update holds at once. + + :param shape: Unused. + :param n_features: Unused. + """ + del shape, n_features + batch_update(self.weights, self.sums, self.counts, self.hx, self.hy) class PerNode: - """Evaluating the neighborhood once per node, which the kernel replaced.""" + """Evaluating the neighborhood once per node, which the contraction replaced.""" params = (SHAPES,) param_names = ("shape",) @@ -124,7 +145,7 @@ def setup(self, shape: tuple[int, int]) -> None: self.nodes = [(x, y) for x in range(shape[0]) for y in range(shape[1])] def time_evaluate_every_node(self, shape: tuple[int, int]) -> None: - """Time the path the kernel replaced, on the same work. + """Time the path the contraction replaced, on the same work. :param shape: Grid shape. """ diff --git a/benchmarks/bench_batch.py b/benchmarks/bench_batch.py deleted file mode 100644 index 1164558..0000000 --- a/benchmarks/bench_batch.py +++ /dev/null @@ -1,160 +0,0 @@ -"""Measure what evaluating the neighborhood once per iteration is worth in batch training. - -Eq. (8) needs ``h_ji`` for every pair of nodes, so a naive loop evaluates the neighborhood once per -node, once per iteration. Because a neighborhood depends only on the offset between two nodes, one -kernel over every offset serves the whole grid and each node's neighborhood is a slice of it. - -Run it directly; it is not part of the test suite, because a timing assertion on shared CI hardware -would be flaky:: - - uv run python benchmarks/bench_batch.py - -Method is the same as ``bench_update.py``, for the same reasons: **interleaved** arms so thermal and -load drift is split evenly rather than attributed to one of them, **medians with an interquartile -range** rather than minima, and **equality asserted first** -- a speed comparison between two -functions that disagree measures nothing. -""" - -from __future__ import annotations - -import functools -import statistics -import timeit -from typing import TYPE_CHECKING - -import numpy as np - -from python_som import SOM -from python_som._core._match import accumulate -from python_som._core._neighborhood import kernel_view, resolve, resolve_kernel -from python_som._core._update import batch_update - -if TYPE_CHECKING: # pragma: no cover - import numpy.typing as npt - -#: Batch iterations per timed run. -ITERATIONS = 12 - -#: Repeats per arm. Odd, so the median is an observation rather than an average of two. -REPEATS = 9 - -#: Grid, sample count and feature count per case. The feature count is varied as well as the grid, -#: because it shifts how much of the work is the contraction rather than the neighborhood, and so -#: changes what there is to win. -CASES = [((20, 20), 200, 4), ((40, 40), 300, 6), ((40, 40), 300, 12), ((60, 60), 400, 8)] - -#: Batch training rejects signed neighborhoods, so only these two can reach this path. -NEIGHBORHOODS = ["gaussian", "bubble"] - -#: Fixed so the reported numbers can be reproduced. -SEED = 20260731 - - -def train( - som: SOM, - data: npt.NDArray[np.floating], - name: str, - *, - use_kernel: bool, -) -> npt.NDArray[np.floating]: - """Run batch training either through the kernel or by evaluating per node. - - Reproduces ``SOM._train_batch`` closely enough to time the difference, rather than calling it, - because the per-node arm no longer exists in the package. - - :param som: A constructed map, used for its shape, radius decay and distance function. - :param data: Training dataset. - :param name: Neighborhood function name. - :param use_kernel: Whether to slice one kernel per iteration or evaluate once per node. - :return: The trained models. - """ - shape = som.get_shape() - weights = som.get_weights().copy() - per_node = resolve(name) - build = resolve_kernel(name) - - for step in range(ITERATIONS): - sigma = som._sigma(step, ITERATIONS) # noqa: SLF001 the decayed radius for this step - sums, counts = accumulate(data, weights, shape, som._distance_function) # noqa: SLF001 - if use_kernel: - kernel = build(shape, sigma, som._cyclic) # noqa: SLF001 - - def neighborhood_of( - node: tuple[int, int], evaluated: npt.NDArray[np.floating] = kernel - ) -> npt.NDArray[np.floating]: - """Slice the kernel for one node.""" - return kernel_view(evaluated, shape, node) - - else: - - def neighborhood_of( - node: tuple[int, int], - radius: float = sigma, - ) -> npt.NDArray[np.floating]: - """Evaluate the neighborhood for one node.""" - return per_node(shape, node, radius, som._cyclic) # noqa: SLF001 - - weights = batch_update(weights, sums, counts, neighborhood_of, shape) - return weights - - -def main() -> None: - """Measure both paths on every case and print the comparison.""" - header = ( - f"{'map':>9} {'samples':>8} {'features':>9} {'h':>12} {'per-node':>11} {'kernel':>11} " - f"{'speedup':>8} {'kernel KB':>10}" - ) - lines = [header, "-" * len(header)] - - for shape, n_samples, n_features in CASES: - rng = np.random.default_rng(SEED) - data = rng.normal(size=(n_samples, n_features)) - - for name in NEIGHBORHOODS: - som = SOM( - x=shape[0], - y=shape[1], - input_len=n_features, - neighborhood_function=name, - neighborhood_radius=3.0, - random_seed=SEED, - ) - som.weight_initialization(mode="random") - - difference = float( - np.abs( - train(som, data, name, use_kernel=True) - - train(som, data, name, use_kernel=False) - ).max() - ) - if difference != 0.0: - message = f"{shape} {name}: the two paths disagree by {difference}" - raise AssertionError(message) - - # functools.partial rather than a lambda: a lambda here would close over the loop - # variables and time whatever they held when it ran, not when it was written. - run_slow = functools.partial(train, som, data, name, use_kernel=False) - run_fast = functools.partial(train, som, data, name, use_kernel=True) - slow: list[float] = [] - fast: list[float] = [] - for _ in range(REPEATS): - slow.append(timeit.timeit(run_slow, number=1)) - fast.append(timeit.timeit(run_fast, number=1)) - - median_slow, median_fast = statistics.median(slow), statistics.median(fast) - kernel_kb = (2 * shape[0] - 1) * (2 * shape[1] - 1) * 8 / 1024 - lines.append( - f"{shape[0]:>4}x{shape[1]:<4} {n_samples:>8} {n_features:>9} {name:>12} " - f"{median_slow * 1e3:>9.1f}ms {median_fast * 1e3:>9.1f}ms " - f"{median_slow / median_fast:>7.2f}x {kernel_kb:>10.0f}" - ) - - lines.append( - f"\nmedian of {REPEATS} interleaved repeats of {ITERATIONS} batch iterations; " - f"both paths verified equal at exactly 0.0 before timing." - ) - print("\n".join(lines)) - - -if __name__ == "__main__": - main() diff --git a/benchmarks/bench_vs_minisom.py b/benchmarks/bench_vs_minisom.py index 9ab2b16..bd5073c 100644 --- a/benchmarks/bench_vs_minisom.py +++ b/benchmarks/bench_vs_minisom.py @@ -5,8 +5,8 @@ uv run python benchmarks/bench_vs_minisom.py -Method is the same as ``bench_update.py`` and ``bench_batch.py``, and the interleaving helper is -imported from the first rather than copied: **interleaved** arms so thermal and load drift is split +Method is the same as ``bench_update.py``, whose interleaving helper is imported rather than +copied: **interleaved** arms so thermal and load drift is split evenly rather than attributed to one of them, **medians with an interquartile range** rather than minima, and **equality asserted first**. @@ -44,8 +44,7 @@ sequential cases that pass is larger than the training itself: 30 steps touch 30 samples, the report touches all 400. The first version of this script did exactly that and reported this package as **7.10x slower** on the largest sequential case. Timing the loop, the same case is **1.08x faster**. -Both numbers were reproducible; only one of them measured training. ``bench_batch.py`` avoids the -same trap by reproducing the loop rather than calling the public method. +Both numbers were reproducible; only one of them measured training. Two tables are printed and the difference between them matters: @@ -368,7 +367,7 @@ def own_initializer() -> list[str]: data = rng.normal(size=(n_samples, n_features)) * 10.0 + 100.0 # functools.partial rather than a closure: a closure over the loop variables would time - # whatever they held when it ran, not when it was written. Same reason bench_batch.py does. + # whatever they held when it ran, not when it was written. ours = functools.partial(seed_and_train_ours, shape, data, n_iteration) theirs = functools.partial(seed_and_train_theirs, shape, data, n_iteration) diff --git a/docs/explanation/batch-vs-stepwise.md b/docs/explanation/batch-vs-stepwise.md index f09ced7..0f653ea 100644 --- a/docs/explanation/batch-vs-stepwise.md +++ b/docs/explanation/batch-vs-stepwise.md @@ -35,13 +35,11 @@ the default is 10 per sample against 1000 for the stepwise modes. from a zeroed array instead destroys them. On a 30×30 map with 20 samples and a small radius, that wiped 282 of 900 models in a single step. -**The per-node sums are contracted with NumPy rather than looped over in Python.** The neighborhood -is evaluated once per node and contracted against the per-node sums and counts. On a 20×20 map with -150 samples this runs about 30× faster than the nested Python loop it replaces, and the two agree -to $10^{-12}$. - -The full $(x, y, x, y)$ tensor would be faster still, but it costs $(xy)^2$ floats, roughly 800 MB -for a 100×100 map, so it is not materialised. +**The sum is contracted, not looped over.** Eq. (8) runs over every pair of nodes, and because the +neighborhood depends only on the offset between two nodes that sum is a convolution. It is evaluated +as two matrix products against per-axis factors, with no loop over nodes and without materialising +the full $(x, y, x, y)$ tensor, which would cost roughly 800 MB for a 100×100 map. See +[How batch training is computed](how-batch-training-is-computed.md). ## Random diff --git a/docs/explanation/comparison-with-som-libraries.md b/docs/explanation/comparison-with-som-libraries.md index aea08dd..df72935 100644 --- a/docs/explanation/comparison-with-som-libraries.md +++ b/docs/explanation/comparison-with-som-libraries.md @@ -3,9 +3,12 @@ Three Python packages implement Kohonen's self-organizing map. This page measures python-som against the other two, and spends most of its length on why that measurement is harder than it looks. -The short version: on batch training python-som is **1.3x to 2.0x faster than SOMPY** and **1.3x to -1.6x slower than MiniSom**, with the gap against MiniSom growing as the map grows. On stepwise -training python-som and MiniSom are within about 10% of each other. SOMPY has no stepwise mode. +The short version, as of 0.7.0: on batch training python-som is **23x to 31x faster than MiniSom** +and **26x to 94x faster than SOMPY**. On stepwise training python-som and MiniSom are within about +10% of each other. SOMPY has no stepwise mode. + +Through 0.6.1 batch training was 1.3x to 1.6x *slower* than MiniSom. +[How batch training is computed](how-batch-training-is-computed.md) covers what changed. ## Why a naive comparison is wrong @@ -74,34 +77,35 @@ With those in place the trained models agree, which is what makes the timings me ## The numbers -Medians of 9 interleaved repeats. python-som 0.6.1, MiniSom 2.3.6 on NumPy 2.5.1, SOMPY @6aca604 on +Medians of 9 interleaved repeats. python-som 0.7.0, MiniSom 2.3.6 on NumPy 2.5.1, SOMPY @6aca604 on NumPy 1.26.4, CPython 3.12.13, Linux, Intel Core Ultra 9 275HX. ### Against MiniSom | map | samples | features | mode | python-som | MiniSom | result | | --- | --- | --- | --- | --- | --- | --- | -| 20x20 | 200 | 4 | batch | 232.5 ms | 239.2 ms | 1.03x faster | -| 40x40 | 300 | 6 | batch | 1016.6 ms | 758.7 ms | **1.34x slower** | -| 60x60 | 400 | 8 | batch | 2885.2 ms | 1763.0 ms | **1.64x slower** | -| 20x20 | 200 | 4 | sequential | 1.2 ms | 1.4 ms | 1.12x faster | -| 40x40 | 300 | 6 | sequential | 2.4 ms | 2.5 ms | 1.06x faster | -| 60x60 | 400 | 8 | sequential | 4.4 ms | 4.7 ms | 1.08x faster | - -**python-som's batch training is slower, and the reason is structural rather than incidental.** -`batch_update` walks every node in Python: on a 60x60 map over 30 iterations that is 108,000 -`einsum` calls. MiniSom's node-side update is a single vectorised divide, and its Python loop runs -over the 400 samples instead. So the ratio tracks nodes against samples, which is why the gap is -absent at 20x20 and 1.64x at 60x60. Vectorising the node loop is the obvious fix, and has not been -done. +| 20x20 | 200 | 4 | batch | 9.8 ms | 227.2 ms | **23.09x faster** | +| 40x40 | 300 | 6 | batch | 22.6 ms | 698.6 ms | **30.93x faster** | +| 60x60 | 400 | 8 | batch | 55.9 ms | 1684.7 ms | **30.15x faster** | +| 20x20 | 200 | 4 | sequential | 1.2 ms | 1.3 ms | 1.11x faster | +| 40x40 | 300 | 6 | sequential | 2.4 ms | 2.6 ms | 1.11x faster | +| 60x60 | 400 | 8 | sequential | 4.1 ms | 4.4 ms | 1.08x faster | + +The batch gap is 0.7.0's doing and comes from two changes: Eq. (8) contracts as two matrix products +instead of one pass per node, and the winner search is one matrix product instead of one norm per +sample. Stepwise is unchanged, and unchanged is the right outcome there: it evaluates one +neighborhood per step, so there was no per-node loop to remove. + +Peak memory is 0.6 MB, 0.9 MB and 1.4 MB against MiniSom's 0.1 MB, 0.5 MB and 1.6 MB. At 100x100 it +is 3.7 MB against 4.6 MB before the change, so the speedup did not cost memory. ### Against SOMPY | map | samples | features | python-som | SOMPY | result | | --- | --- | --- | --- | --- | --- | -| 20x20 | 200 | 4 | 145.4 ms | 188.3 ms | 1.30x faster | -| 40x40 | 300 | 6 | 727.7 ms | 1264.6 ms | 1.74x faster | -| 60x60 | 400 | 8 | 1906.4 ms | 3796.6 ms | 1.99x faster | +| 20x20 | 200 | 4 | 6.7 ms | 176.1 ms | **26.20x faster** | +| 40x40 | 300 | 6 | 17.1 ms | 1201.3 ms | **70.31x faster** | +| 60x60 | 400 | 8 | 38.6 ms | 3642.1 ms | **94.36x faster** | Batch only, since SOMPY implements nothing else: `SOMFactory.build` accepts `training='seq'` and ignores it. @@ -121,7 +125,7 @@ scaled by 10 and offset 100 from the origin: | --- | --- | --- | --- | | 20x20 | 1.6797 | 1.1157 | MiniSom | | 40x40 | 0.2855 | 0.6265 | python-som | -| 60x60 | 0.0937 | 0.5009 | python-som | +| 60x60 | 0.0873 | 0.5009 | python-som | The three initializers place models on the plane of the first two principal components and differ in how far apart. python-som uses $\bar{x} + c_1\sqrt{\lambda_1}v_1 + c_2\sqrt{\lambda_2}v_2$, scaling @@ -182,7 +186,8 @@ because shared runners cannot measure anything: an earlier version of this compa The peers are pinned, and MiniSom's development branch has already moved past its release with a numba-compiled batch path that would change these numbers substantially. -Both peers are worth using. MiniSom is faster at batch training on larger maps and has hexagonal -topologies, which python-som does not. SOMPY has clustering and visualization helpers built in. What -python-som offers against them is NumPy as its only runtime dependency, a scikit-learn estimator -interface, pickle-free artifacts with provenance, and type information. +Both peers are worth using. MiniSom has hexagonal topologies, a triangular neighborhood and several +distance metrics, none of which python-som has. SOMPY has clustering and visualization helpers built +in. What python-som offers against them is speed on batch training, NumPy as its only runtime +dependency, a scikit-learn estimator interface, pickle-free artifacts with provenance, and type +information. diff --git a/docs/explanation/how-batch-training-is-computed.md b/docs/explanation/how-batch-training-is-computed.md new file mode 100644 index 0000000..48c1119 --- /dev/null +++ b/docs/explanation/how-batch-training-is-computed.md @@ -0,0 +1,144 @@ +# How batch training is computed + +Kohonen's Eq. (8) says what a batch update is. It does not say how to evaluate it, and the +difference is a factor of thirty. + +$$m_i^* = \frac{\sum_j n_j h_{ji} \bar{x}_{m,j}}{\sum_j n_j h_{ji}}$$ + +Read literally, that is a sum over every pair of nodes, evaluated once per node. On a 60x60 map +over 30 iterations, the literal reading is 108,000 evaluations of the neighborhood. This package +does it in two matrix products. + +## The sum is a convolution + +$h_{ji}$ depends only on the offset between nodes $j$ and $i$, never on where either sits. That is +what makes the map translation-invariant, and it means the numerator is a convolution of the +per-node sums with the neighborhood, and the denominator a convolution of the per-node counts. + +A convolution can be evaluated many ways. The one that wins here depends on a second property. + +## Both batch neighborhoods are separable + +Batch training admits the gaussian and the bubble, and each factors into a product of per-axis +terms: + +$$e^{-(dx^2 + dy^2) / 2\sigma^2} = e^{-dx^2 / 2\sigma^2} \cdot e^{-dy^2 / 2\sigma^2}$$ + +$$\max(|dx|, |dy|) \le r \iff (|dx| \le r) \land (|dy| \le r)$$ + +The first is a property of the exponential. The second is a property of the Chebyshev metric, which +is the metric this package's bubble uses; a Euclidean disc would not factor. Neither is a property +of neighborhood functions in general, and the mexican hat has no such factorisation, which is one +of two reasons batch training rejects it. + +Given the factors as matrices $H^x_{ac} = f(a-c)$ and $H^y_{bd} = g(b-d)$, the whole update is: + +```python +numerator = np.einsum("ac,bd,cdf->abf", hx, hy, sums, optimize=True) +denominator = np.einsum("ac,bd,cd->ab", hx, hy, counts, optimize=True) +``` + +$H^x$ is $X \times X$ and $H^y$ is $Y \times Y$, so the memory is $X^2 + Y^2$ floats: 58 KB on a +60x60 map, against the 104 MB a full node-by-node matrix would need and the 800 MB it would need at +100x100. + +Measured against evaluating the neighborhood per node: + +| map | per node | axis matrices | | +| --- | --- | --- | --- | +| 20x20 | 1.79 ms | 0.054 ms | 33x | +| 40x40 | 14.64 ms | 0.093 ms | 158x | +| 60x60 | 65.11 ms | 0.135 ms | 482x | + +## This is not the separability mistake + +The distinction matters, because the two look identical from a distance and this package shipped +the wrong one once. + +An **axis profile** is a way of evaluating a neighborhood that is already defined as a function of +$\mathrm{sqdist}$. The definition does not change; only the order of the arithmetic does, and a test +asserts the outer product of the two factors equals the isotropic function node by node. + +A **separably defined neighborhood** is a different function. Building a mexican hat as an outer +product of two one-dimensional Ricker wavelets gives $+0.165$ on the diagonal at $2\sigma$ where the +correct value is $-0.055$: an excitatory lobe exactly where the function must inhibit. That was a +real defect here, and [Why isotropy matters](why-isotropy-matters.md) covers it. + +The guard is a registry. A neighborhood has an axis profile only where the factorisation is an +identity, and a test asserts the registry holds exactly the unsigned neighborhoods. A future +neighborhood that is unsigned but not separable fails that test rather than being approximated. + +## Finding every winner at once + +The update is now a small part of the cost. The larger part is Eq. (4), the search for each sample's +best-matching model: + +$$c = \arg\min_i \lVert x - m_i \rVert$$ + +Expanding the norm gives $\lVert x \rVert^2 - 2\,x \cdot w + \lVert w \rVert^2$, and the first term +is the same for every model, so it cannot change which one wins. What remains is a matrix product +against all the models at once, plus a per-node constant. + +**This is not Kohonen's dot-product map.** Section 4.5 defines a genuinely different algorithm, +$c = \arg\max_i \mathrm{dot}(x, m_i)$, which requires the models to be renormalized to constant +length after every cycle and picks a different node when they are not. The expansion above is exact +for the Euclidean distance and needs no normalization. + +### The expansion cancels, and the fix is one line + +$\lVert w \rVert^2$ grows with the square of the data's distance from the origin, while the +differences between models do not. With models offset by $10^9$, that term is around $10^{18}$ and +the subtraction loses every significant digit: + +| offset | samples given the wrong node | +| --- | --- | +| origin, 1e3, 1e6 | 0 of 500 | +| **1e9** | **500 of 500** | +| **1e12** | **500 of 500** | + +Subtracting a common shift from both sides is exact in $\lVert x - w \rVert$, costs 1%, and removes +it at every offset tested. Data far from the origin is not exotic: timestamps, easting and northing +coordinates and absolute sensor readings all look like this. It is the same failure mode that +[linear initialization](why-linear-initialization-is-an-svd.md) had before 0.4.0. + +A custom `distance_function` keeps the exact per-sample loop, because the expansion is an identity +for the Euclidean norm and nothing else. + +### Small blocks beat large ones + +The search runs in blocks so the score matrix never grows with the dataset. The block size is +tuned rather than chosen, on a 60x60 map with 2000 samples: + +| budget | time | peak | +| --- | --- | --- | +| **512 KB** | **7.62 ms** | **1.07 MB** | +| 2 MB | 7.50 ms | 2.57 MB | +| 8 MB | 11.14 ms | 8.56 MB | + +A block that fits in cache is read back by `argmin` for free. One that does not is read back from +memory, which is why the largest budget is both the slowest and the heaviest. + +## What Kohonen says about all this + +Reorganising the arithmetic is not a departure from the paper. Section 4.4 derives Eq. (8) from +Eq. (7) on exactly these grounds, that "the same addends occur a great number of times", and +Section 5.2 notes that Eq. (8) "allows for a very efficient implementation" and that "the winner +search can be partly parallelized by dividing the data". + +One requirement does constrain the implementation. Section 4.4 closes: the old values "are replaced +by the respective means, **in one concurrent computing operation over all nodes of the grid**". +Every node must be computed from the models as they stood at the start of the iteration. The +contraction satisfies this structurally, since there is no loop to get wrong, and a test asserts it +directly. + +Two further optimizations the paper suggests are **not** implemented here. Section 5.2 proposes +confining the winner search to the neighborhood of the previous winner, which is an approximation +that can miss a better match, and reducing the models to eight-bit precision, which changes results +by far more than round-off. Both are recorded rather than adopted. + +## What this cost + +Trained weights differ from 0.6.1 by about $10^{-15}$ relative. The contraction sums the same terms +in a different order, so the results are not bit-identical to earlier versions, and +[Reproduce a result](../how-to/reproduce-a-result.md) says which version to pin to reproduce an +older figure exactly. diff --git a/docs/explanation/why-isotropy-matters.md b/docs/explanation/why-isotropy-matters.md index 9ebc9d5..e6aa0dc 100644 --- a/docs/explanation/why-isotropy-matters.md +++ b/docs/explanation/why-isotropy-matters.md @@ -61,8 +61,14 @@ requires. ## How the package enforces it Every neighborhood function is built from `squared_grid_distance`, which reduces the two offsets to -one number before any profile is applied. The three shipped functions share a single implementation -of each formula, so the per-node form and the batch kernel cannot drift apart. +one number before any profile is applied. + +Batch training evaluates the same functions by a different route, contracting per-axis factors +instead of calling them once per node. That is a contraction strategy rather than a second +definition, and it is available only for the two neighborhoods where the factorisation is an +identity. A test asserts the factors multiply back to the isotropic function node by node, so the +two cannot drift apart, and the mexican hat has no factor at all. See +[How batch training is computed](how-batch-training-is-computed.md). The test suite asserts the property directly rather than checking golden values: equal grid distance must give equal $h$. That assertion fails against the separable construction and passes against the @@ -76,9 +82,9 @@ rather than a disc. Kohonen's phrasing ("up to a certain radius from the winner" so the two sources genuinely differ; this package follows Vrieze and says so rather than quietly picking one. -The consequence is worth stating because it is easy to assume otherwise: on a large enough grid, -nodes at equal Euclidean distance can fall on opposite sides of the boundary. The smallest case is a -radius of $\sqrt{50}$, where $(5, 5)$ lies inside a $\sigma = 5$ square and $(7, 1)$ lies outside. +A Chebyshev ball is not isotropic under the Euclidean metric. On a large enough grid, nodes at equal +Euclidean distance can fall on opposite sides of the boundary: at a radius of $\sqrt{50}$, $(5, 5)$ +lies inside a $\sigma = 5$ square and $(7, 1)$ lies outside. ## Further reading diff --git a/docs/explanation/why-linear-initialization-is-an-svd.md b/docs/explanation/why-linear-initialization-is-an-svd.md new file mode 100644 index 0000000..35bd541 --- /dev/null +++ b/docs/explanation/why-linear-initialization-is-an-svd.md @@ -0,0 +1,87 @@ +# Why linear initialization is an SVD + +Kohonen recommends starting the models on the plane of the data's two largest principal components +rather than at random, because "much faster convergence follows" (Section 4.3). Computing those +components is the only linear algebra this package needs. + +Through 0.3.0 it was `sklearn.decomposition.PCA`. Since 0.4.0 it is about twenty lines of +`np.linalg.svd`, which removed a dependency and also fixed an accuracy defect. + +## Two ways to find the same components + +For centred data $X$, the principal components are the eigenvectors of the covariance matrix +$X^\top X / (n-1)$. There are two ways to get them. + +**Eigendecompose the covariance matrix.** Form $X^\top X$, then decompose it. Cheap when there are +far more samples than features. + +**Decompose the data directly.** For $X = U S V^\top$, the rows of $V^\top$ are the components and +$S^2/(n-1)$ the variance along each. No covariance matrix is ever formed. + +They agree in exact arithmetic. They do not agree in floating point, because forming $X^\top X$ +**squares the condition number**. Every digit of precision in the data becomes half a digit in the +result, and when the mean is large relative to the spread there are not many digits to start with. + +## The defect this exposed + +Linear initialization fits its PCA on **raw** data by design, so the models live in the same space +as the inputs they will be compared against. Since scikit-learn 1.5 the default solver picks +`covariance_eigh` when samples comfortably outnumber features, which is exactly the squaring path. + +On `(150, 4)` data offset by $10^7$, the second explained variance was wrong by **5.8%**, and the +models it produced differed from the correct ones by 2.43 against a total model spread of 2.0. The +error was larger than the structure being initialized. + +Measured against a reference centred in `longdouble` before decomposing: + +| solver | relative error | +| --- | --- | +| scikit-learn `auto` (covariance path) | 1.4e-06 to 5.5e-06 | +| scikit-learn `svd_solver="full"` | ~1e-15 | +| this package | ~1e-15 | + +Data offset far from the origin is not a corner case. Timestamps, easting and northing coordinates +and absolute sensor readings all look like this, and none of them announce themselves. + +The same failure mode appears again in the best-matching-unit search, where expanding +$\lVert x - w \rVert^2$ squares the magnitudes in the same way. +[How batch training is computed](how-batch-training-is-computed.md) covers it and the one-line fix. + +## Why the reimplementation is trusted + +Replacing a widely-used library's numerics with twenty lines of your own is the change in this +package a reviewer should be least willing to take on faith, so it is not asked for on faith. + +scikit-learn remains a **test** dependency, and `tests/test_linalg_matches_sklearn.py` re-derives +every fit both ways on every CI run and compares them. The claim under test is not "close enough" +but "the same numbers": the tolerances are at the scale of double-precision round-off. + +The comparison is against `svd_solver="full"` rather than the default, because the default is the +inaccurate path and comparing against it would fail a correct implementation. A second check uses a +`longdouble` reference, which depends on no library's solver choice and would survive scikit-learn +changing its defaults again. + +## Two details that are easy to get wrong + +**The sign convention is v-based.** An SVD fixes each component only up to sign, so a convention is +needed for a fit to be reproducible. scikit-learn's PCA calls `svd_flip` with +`u_based_decision=False`, orienting each component so its largest-magnitude loading is positive. +That is the less common of the two settings in that helper. Taking the default would still give a +valid PCA, but a different one, and linear initialization would lay its models out reversed along +that axis. No test of orthonormality or explained variance would notice; only comparing signs does. + +**A near-constant column is scaled by 1.** Dividing a column by its own standard deviation is the +obvious z-score and the wrong one when that deviation is zero. The guard is not `variance == 0` +either: a column built by arithmetic that should cancel exactly can retain a variance of about +$10^{-30}$, which passes an equality test and then divides the column by roughly $10^{-15}$. The +bound used is scikit-learn's, from Chan, Golub and LeVeque. + +## What this costs + +The reimplementation removed 264 MB of required install, 79% of the payload, taking python-som from +10 packages to 1. That was the reason for doing it. The accuracy improvement was a side effect, and +in retrospect the more valuable half. + +Because 0.4.0 changed what linear initialization produces for data far from the origin, results are +not comparable with 0.3.0 for those datasets. +[Reproduce a result](../how-to/reproduce-a-result.md) has the version-pinning details. diff --git a/docs/how-to/reproduce-a-result.md b/docs/how-to/reproduce-a-result.md index 73438d7..9597648 100644 --- a/docs/how-to/reproduce-a-result.md +++ b/docs/how-to/reproduce-a-result.md @@ -25,16 +25,19 @@ Numerical results are allowed to change between minor versions before 1.0.0, and versions after it. A seed alone does not pin a result across an upgrade. ``` -python-som==0.4.0 +python-som==0.7.0 ``` -Two specific breaks worth knowing about, if you are reproducing an older figure: +Three releases change results. If you are reproducing an older figure: - **0.3.0** replaced the global RNG with a per-instance generator, so `random_seed=42` gives a different map from 0.2.0 and earlier. Pin `python-som==0.2.0` to reproduce those. - **0.4.0** fixed an accuracy defect in linear initialization for data far from the origin. Near the origin the difference is floating-point noise; far from it, it is large, and 0.4.0 is the correct one. +- **0.7.0** made batch training 20x to 40x faster by reorganising the same arithmetic, which sums in + a different order. Weights differ from 0.6.1 by about 1e-15 relative: far below anything a result + depends on, and still not bit-identical. Pin `python-som==0.6.1` if you need an exact match. ## Record what you ran @@ -48,7 +51,7 @@ print(som.last_report) ``` TrainingReport(mode='batch', n_iteration=100, n_samples=150, random_seed=42, final_learning_rate=None, final_neighborhood_radius=0.5, - quantization_error=0.3142, python_som_version='0.4.0', + quantization_error=0.3142, python_som_version='0.7.0', numpy_version='2.5.1', wall_time_seconds=0.42) ``` diff --git a/docs/how-to/save-and-load-a-map.md b/docs/how-to/save-and-load-a-map.md index 18d645d..1d0e7e8 100644 --- a/docs/how-to/save-and-load-a-map.md +++ b/docs/how-to/save-and-load-a-map.md @@ -32,7 +32,7 @@ with np.load("iris-map.npz", allow_pickle=False) as archive: ```json { "format_version": 1, - "python_som_version": "0.4.0", + "python_som_version": "0.7.0", "numpy_version": "2.5.1", "config": { "shape": [10, 10], "input_len": 4, @@ -48,7 +48,7 @@ with np.load("iris-map.npz", allow_pickle=False) as archive: "mode": "batch", "n_iteration": 100, "n_samples": 150, "random_seed": 42, "final_learning_rate": null, "final_neighborhood_radius": 0.5, "quantization_error": 0.3142, - "python_som_version": "0.4.0", "numpy_version": "2.5.1", "wall_time_seconds": 0.42 + "python_som_version": "0.7.0", "numpy_version": "2.5.1", "wall_time_seconds": 0.42 } } ``` diff --git a/docs/how-to/use-a-custom-strategy.md b/docs/how-to/use-a-custom-strategy.md index f2f21ce..a3e4ee9 100644 --- a/docs/how-to/use-a-custom-strategy.md +++ b/docs/how-to/use-a-custom-strategy.md @@ -30,9 +30,8 @@ between two nodes, not of the two axis offsets separately. See ## What a custom strategy costs you -One thing, and it is worth knowing before you commit: a callable cannot be written to a file without -`pickle`, so `save_npz` records only its **name**. A map trained with your own function will not -reload on its own. +A callable cannot be written to a file without `pickle`, so `save_npz` records only its **name**. A +map trained with your own function will not reload on its own. ```python python_som.SOM.load_npz("custom.npz") diff --git a/docs/how-to/use-with-scikit-learn.md b/docs/how-to/use-with-scikit-learn.md index 903fda4..6e7febc 100644 --- a/docs/how-to/use-with-scikit-learn.md +++ b/docs/how-to/use-with-scikit-learn.md @@ -39,7 +39,7 @@ rows, columns = np.unravel_index(som.predict(X), som.get_shape()) `winner(x)` still returns `(row, column)` for a single sample. -## Two differences from scikit-learn worth knowing +## Two differences from scikit-learn **`fit` continues rather than resetting.** scikit-learn estimators conventionally discard their fitted state on a second `fit`. This one does not: a SOM's models *are* its state, and `train` has diff --git a/docs/reference/api.md b/docs/reference/api.md index 93ca9c0..296dc05 100644 --- a/docs/reference/api.md +++ b/docs/reference/api.md @@ -14,12 +14,11 @@ - bubble - mexican_hat - axis_offsets - - offset_span + - axis_matrix + - gaussian_axis_profile + - bubble_axis_profile + - resolve_axis_profile - squared_grid_distance - - gaussian_kernel - - bubble_kernel - - mexican_hat_kernel - - kernel_view - resolve ## Decay functions diff --git a/docs/reference/neighborhood-functions.md b/docs/reference/neighborhood-functions.md index 549df84..a688786 100644 --- a/docs/reference/neighborhood-functions.md +++ b/docs/reference/neighborhood-functions.md @@ -1,8 +1,8 @@ # Neighborhood functions -The neighborhood function decides how a winner's correction spreads to the rest of the grid. It is -the part of the SOM that turns vector quantization into a *topology-preserving* map, and the part -where a plausible-looking implementation can be quietly wrong, so it is worth setting out in full. +The neighborhood function decides how a winner's correction spreads to the rest of the grid. This +page is the formulas and constants; [Why isotropy matters](../explanation/why-isotropy-matters.md) +covers why they take the form they do. ## The rule that governs all of them @@ -100,16 +100,10 @@ sometimes even better"* than a distance-dependent one. !!! note "The bubble uses the Chebyshev metric, so the region is a square" - A node is included when $\max(|dx|, |dy|) \le \rho$, which makes the region a square rather than - a disc. That follows Vrieze's appendix pseudo-code, which computes - `b = MAX(ABS(i - w_i), ABS(j - w_j))`, although Kohonen's phrase "up to a certain radius from the - winner" reads as Euclidean. The two sources genuinely differ, and this library keeps Vrieze's - reading so that existing results stay reproducible. - - A consequence worth knowing, because it is easy to assume otherwise: **a Chebyshev ball is not - isotropic under the Euclidean metric**. Two nodes the same Euclidean distance from the winner can - fall on opposite sides of the boundary. The smallest case is $r = \sqrt{50}$, where $(5, 5)$ is - inside a $\sigma = 5$ square and $(7, 1)$ is outside. + A node is included when $\max(|dx|, |dy|) \le \rho$, following Vrieze's appendix pseudo-code, + `b = MAX(ABS(i - w_i), ABS(j - w_j))`. Kohonen's phrase "up to a certain radius from the winner" + reads as Euclidean, so the two sources differ; the consequences are in + [Why isotropy matters](../explanation/why-isotropy-matters.md). Unlike the other two, a radius of zero is allowed here: it selects the winner alone, which is well defined for an indicator function, where for the gaussian and the mexican hat it would be a division diff --git a/mkdocs.yml b/mkdocs.yml index a043061..d52e3c9 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -43,6 +43,8 @@ nav: - Changelog: reference/changelog.md - Explanation: - Why isotropy matters: explanation/why-isotropy-matters.md + - How batch training is computed: explanation/how-batch-training-is-computed.md + - Why linear initialization is an SVD: explanation/why-linear-initialization-is-an-svd.md - Batch vs stepwise: explanation/batch-vs-stepwise.md - Artifact safety: explanation/artifact-safety.md - Comparison with MiniSom and SOMPY: explanation/comparison-with-som-libraries.md diff --git a/pyproject.toml b/pyproject.toml index 8bf3e98..a8cd067 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "hatchling.build" [project] name = "python-som" -version = "0.6.1" +version = "0.7.0" authors = [{ name = "André Moreira Souza", email = "msouza.andre@hotmail.com" }] description = "Python implementation of the Self-Organizing Map" readme = "README.md" @@ -26,13 +26,9 @@ classifiers = [ "Operating System :: OS Independent", "Typing :: Typed", ] -# NumPy is the only thing this package needs at runtime. pandas and scikit-learn were dropped in -# 0.4.0: pandas because `np.asarray` already converts a DataFrame through the `__array__` protocol, -# and scikit-learn because its PCA and StandardScaler are twenty lines of `np.linalg.svd`. Together -# together they pulled in 8 further packages to reach that. Measured on Linux/CPython 3.12, the -# installed footprint drops from 333 MB across 10 packages to 69 MB across 1: a 264 MB reduction, -# 79% of the payload. Both remain in -# `dev`, where the differential tests re-check the replacements against them on every CI run. +# NumPy is the only runtime dependency. pandas and scikit-learn were dropped in 0.4.0, taking the +# install from 333 MB across 10 packages to 69 MB across 1. Both stay in `dev`, where differential +# tests re-check the replacements against them. dependencies = [ "numpy>=1.24", ] @@ -41,25 +37,18 @@ dependencies = [ cli = ["tqdm>=4.66"] dev = [ "hypothesis==6.163.0", - # Not a runtime dependency. tests/test_minisom_agreement.py checks our neighborhood functions, - # decay and training loops against MiniSom's, the closest comparable implementation, the same way - # tests/test_linalg_matches_sklearn.py checks our PCA against scikit-learn's. Two independent - # implementations of Kohonen's equations agreeing is worth more than either one's own tests. - # Cheap to depend on: MIT, a single minisom.py, and it declares no dependencies of its own. + # Test-only. tests/test_minisom_agreement.py checks this package's equations against MiniSom's. + # MIT, a single file, and it declares no dependencies of its own. "minisom==2.3.6", "mypy==2.3.0", - # Not runtime dependencies. The suite exercises the DataFrame path through the port, and - # tests/test_linalg_matches_sklearn.py re-derives every PCA both ways and compares. Pinned like - # the rest of the tooling: scikit-learn changed its default PCA solver in 1.5, and a floating pin - # would surface that as an unrelated PR's CI failure instead of a Dependabot diff to read. - # Split by Python version because the current majors dropped 3.10, which this package still - # supports: same reason as tomli below, same pattern. + # Test-only: the DataFrame path through the port, and the differential PCA test. Pinned because + # scikit-learn changed its default PCA solver in 1.5, and a floating pin would surface that as + # an unrelated PR's CI failure. Split by Python version, since the current majors dropped 3.10. "pandas==3.0.5; python_version >= '3.11'", "pandas==2.3.3; python_version < '3.11'", "scikit-learn==1.9.0; python_version >= '3.11'", "scikit-learn==1.7.2; python_version < '3.11'", - # pandas-stubs 3.x needs Python >=3.11; it is dev-only, so gate it rather than - # raising the library's own floor. mypy runs on 3.12 in CI, where it is present. + # pandas-stubs 3.x needs Python >=3.11; gated rather than raising the library's floor. "pandas-stubs==3.0.3.260530; python_version >= '3.11'", "pre-commit==4.6.1", "pytest==9.1.1", @@ -72,48 +61,30 @@ dev = [ ] docs = [ "mkdocs-material==9.7.7", - # Pinned to 1.2.2 deliberately, not to the latest 1.2.3. The canonical repository, - # github.com/mkdocs/mkdocs-redirects (245 stars, active, not archived), has tags up to v1.2.2 and - # no v1.2.3. PyPI's 1.2.3 was published on 2026-03-28 declaring its source as - # github.com/ProperDocs/properdocs-redirects, a repository created 2026-03-14 with 3 stars -- a - # version upstream never tagged, released from a repository that appeared two weeks earlier while - # the original stayed active. 1.2.2's declared source is the canonical repository and matches its - # own v1.2.2 tag, so that is the last release whose provenance can be checked. - # - # Do not bump this without re-checking who publishes it. + # SUPPLY CHAIN: pinned to 1.2.2, not the newer 1.2.3. Upstream (github.com/mkdocs/mkdocs- + # redirects) has no v1.2.3 tag; PyPI's 1.2.3 declares its source as a different repository + # created two weeks before it was published. 1.2.2 is the last release whose provenance checks + # out. Do not bump without re-checking who publishes it. "mkdocs-redirects==1.2.2", "mkdocstrings-python==2.0.5", ] -# Static analysis used during review, kept out of `dev` so the CI jobs that never run it do not -# install it. Pinned exactly, like the rest of the tooling. Nothing here gates a merge: ruff and mypy -# are the enforced checks, and these are for the deeper passes a human asks for. +# Review tools, kept out of `dev` so CI jobs that never run them do not install them. Nothing here +# gates a merge; ruff and mypy are the enforced checks. # -# `osv-scanner` is deliberately absent. The real tool is a Go binary from Google; the PyPI package of -# that name was registered on 2026-07-16 with one release, no author, no home page, and a summary of -# "Reserved name placeholder. No functionality." Install the Go binary if you want it. +# SUPPLY CHAIN, three tools deliberately absent: +# osv-scanner the PyPI package of that name is a placeholder with no functionality; the real +# tool is a Go binary from Google. +# guarddog genuine (DataDog) but first released three months after this project, so it fails +# the rule that a dependency should predate what it is added to. +# semgrep pins mcp==1.23.3 exactly, which carries PYSEC-2026-3481/3482/3483 and so cannot be +# remediated from here. Revisit if it relaxes the pin. +# Benchmarking, kept out of `dev`: asv pulls about ten transitive packages and no gating job runs +# it. minisom is repeated so `--extra bench` alone runs every script in benchmarks/. # -# `guarddog` is also absent: genuine (DataDog) but first released 2022-11-28, three months *after* -# this project, so it fails the rule that a dependency should predate what it is added to. Add it -# deliberately if the supply-chain checks are wanted, not as part of a batch. -# -# `semgrep` was added and then removed. It pins `mcp==1.23.3` exactly, and that version carries -# PYSEC-2026-3481/3482/3483 (fixed in mcp 1.27.2/1.28.1) -- so the exact pin makes the advisories -# unremediable from here, and 1.172.0 is already the latest semgrep. Against that, bandit reports -# nothing in `src/`, and this is a pure-numpy library with no network, auth, crypto, SQL or -# templating, which is most of what semgrep's rulesets look for. Low coverage gained for three CVEs -# that cannot be patched. Revisit if semgrep relaxes the pin. -# Benchmarking, kept out of `dev` for the same reason `analysis` is: `asv` pulls about ten -# transitive packages (asv-runner, json5, build, tabulate, virtualenv, packaging, -# importlib-metadata, pyyaml, pympler) and no gating CI job runs it. `minisom` is repeated here so -# that `--extra bench` alone is enough to run every script in benchmarks/. -# -# SOMPY is deliberately absent and cannot be added. It uses `np.Inf`, removed in NumPy 2.0, at -# class-definition time, so it cannot be imported in any environment this package supports. -# `benchmarks/bench_vs_sompy.py` drives it through a separate interpreter instead; that environment -# is built by hand and the command is in the script's docstring. +# SOMPY cannot be added: it uses `np.Inf`, removed in NumPy 2.0, at class-definition time. +# benchmarks/bench_vs_sompy.py drives it through a separate interpreter built by hand. bench = [ - # airspeed velocity, from the airspeed-velocity organisation: 22 releases since 2015-05-01, - # and the tool numpy, scipy, pandas and scikit-learn all use to track their own performance. + # What numpy, scipy, pandas and scikit-learn all use to track their own performance. "asv==0.6.6", "minisom==2.3.6", ] @@ -126,10 +97,8 @@ analysis = [ # reports can be reproduced from the command line. "pylint==4.0.6", ] -# Restores the install set that 0.3.0 pulled in by default, for anyone who was relying on it. -# scikit-learn integration. The core never imports it; only python_som.sklearn does, and that module -# exists because these paths need `__sklearn_tags__`, which in practice means inheriting BaseEstimator. -# A lower bound rather than a pin: this is a user-facing install set, not tooling. +# scikit-learn integration. Only python_som.sklearn imports it. Lower bounds rather than pins: +# these are user-facing install sets, not tooling. `examples` restores what 0.3.0 installed. sklearn = [ "scikit-learn>=1.4", ] @@ -187,9 +156,8 @@ ignore = [ ] [tool.ruff.lint.per-file-ignores] -# Tests must be able to build a DataFrame to exercise the port, and to compare our PCA against -# sklearn's differentially. The ban protects the core, not the suite that checks it. -# The scikit-learn adapter is the module that exists to import scikit-learn. +# The ban protects the core, not the suite that checks it, nor the adapter that exists to import +# scikit-learn. "src/python_som/sklearn.py" = ["TID251"] "tests/test_sklearn_adapter.py" = ["TID251"] "tests/*" = [ @@ -205,14 +173,12 @@ ignore = [ "PLR2004", ] "examples/*" = ["INP001", "T201", "D100"] -# Scripts a human runs, not importable modules: a table on stdout is the deliverable, and neither -# directory should become a package just to satisfy the implicit-namespace rule. +# Scripts a human runs: a table on stdout is the deliverable, and they are not importable modules. "benchmarks/*" = ["INP001", "T201"] [tool.ruff.lint.flake8-tidy-imports.banned-api] -# The package is numpy-only at runtime, and as of 0.4.0 this ban has no per-file exemption anywhere -# under src/: the two modules that used to need one no longer import either library. Only tests are -# exempt, so the replacements can be compared against the originals. +# The package is numpy-only at runtime. No module under src/ is exempt; only tests are, so the +# replacements can be compared against the originals. "pandas".msg = "The core is numpy-only. Convert at the boundary in python_som/_convert.py." "sklearn".msg = "The core is numpy-only. Linear algebra belongs in python_som/_core/_linalg.py." @@ -223,26 +189,32 @@ convention = "pep257" known-first-party = ["python_som"] [tool.mypy] -# Deliberately not pinned to python_version = "3.10". NumPy's bundled stubs use 3.12 syntax, and -# `tomllib` is 3.11+, so pinning to the floor makes mypy fail on the dependencies rather than on -# our code. Real 3.10 compatibility is covered by the CI test matrix instead. +# Not pinned to 3.10: NumPy's stubs use 3.12 syntax, so pinning to the floor makes mypy fail on the +# dependencies rather than on our code. Real 3.10 support is covered by the CI matrix. strict = true files = ["src", "tests"] warn_unreachable = true enable_error_code = ["ignore-without-code", "redundant-expr", "truthy-bool"] [[tool.mypy.overrides]] -module = ["sklearn.*"] +module = ["sklearn.*", "numba.*"] ignore_missing_imports = true -# scikit-learn ships no py.typed, so mypy sees BaseEstimator as Any and --strict refuses to subclass -# it. Relaxed for the adapter module alone, which is the only place that inherits from scikit-learn; -# everything else stays under full strictness. The alternative, a blanket type: ignore on the class, -# would hide any *other* subclassing mistake in the same file. +# scikit-learn ships no py.typed, so --strict refuses to subclass BaseEstimator. Relaxed for the +# adapter alone; a blanket type: ignore would hide any other mistake in the same file. [[tool.mypy.overrides]] module = ["python_som.sklearn"] disallow_subclassing_any = false +# numba ships no py.typed, so --strict rejects the jit decorator and `prange` as untyped. Relaxed +# for the accelerator alone; its contract is the BmuKernel protocol, checked by a differential test. +# These only fire when numba is installed, which is the accelerated-path CI job. +[[tool.mypy.overrides]] +module = ["python_som._accelerate"] +disallow_untyped_decorators = false +disallow_untyped_calls = false +disable_error_code = ["attr-defined"] + [tool.pytest.ini_options] minversion = "8.0" testpaths = ["tests"] @@ -257,8 +229,7 @@ markers = [ filterwarnings = ["error"] [tool.coverage.run] -# source_pkgs, not source: with a src layout the package under test is the installed one, and -# coverage has to resolve it by import name rather than by directory. +# source_pkgs, not source: with a src layout the package under test is the installed one. source_pkgs = ["python_som"] branch = true @@ -266,23 +237,19 @@ branch = true source = ["src/python_som", "*/site-packages/python_som"] [tool.coverage.report] -# The suite is at 100%. A lower gate would permit a silent regression, and anything genuinely -# unreachable should carry an explicit `# pragma: no cover` a reviewer can question. +# Anything genuinely unreachable carries an explicit `# pragma: no cover` a reviewer can question. fail_under = 100 show_missing = true exclude_lines = [ "pragma: no cover", "if TYPE_CHECKING:", "raise NotImplementedError", - # A bare `...` is a Protocol method body: a declaration of shape that is never executed, since - # the protocols are structural and nothing inherits from them. Excluded as a rule rather than - # with four identical pragmas, and narrow enough to stay honest -- it matches only a line whose - # entire content is an ellipsis, which in this package occurs solely in _core/_protocols.py. + # A bare `...` is a Protocol method body, never executed. Matches only a line whose entire + # content is an ellipsis, which occurs solely in _core/_protocols.py. "^\\s*\\.\\.\\.$", ] [tool.bandit] -# `assert` is how a test states its expectation, so B101 fires on every one of them and says nothing. -# Excluded by directory rather than by skipping the check, so B101 still applies to `src/`, where an -# assert *would* be wrong: asserts vanish under `python -O`, so validation must raise instead. +# B101 fires on every test assert. Excluded by directory rather than skipped, so it still applies +# to `src/`, where an assert would be wrong: they vanish under `python -O`. exclude_dirs = ["tests", ".venv", "site", "dist"] diff --git a/src/python_som/__init__.py b/src/python_som/__init__.py index 1bd3f1c..7e0c852 100644 --- a/src/python_som/__init__.py +++ b/src/python_som/__init__.py @@ -9,10 +9,8 @@ >>> som.weight_initialization(mode="linear", data=data) >>> error = som.train(data, n_iteration=100, mode="batch") -Internally the package is a pure functional core with a thin shell around it: -:mod:`python_som._core` holds every numeric decision as functions over NumPy arrays, and imports -nothing but NumPy. :mod:`python_som._convert` adapts pandas and anything else array-like at the -boundary, and :mod:`python_som._som` holds the state and the training loops. +Internally, :mod:`python_som._core` holds every numeric decision as pure functions over NumPy +arrays and imports nothing else; a thin shell around it handles conversion, state and I/O. Reference: Teuvo Kohonen, Essentials of the self-organizing map, Neural Networks 37 (2013) 52-65, diff --git a/src/python_som/_accelerate.py b/src/python_som/_accelerate.py new file mode 100644 index 0000000..1a7ca55 --- /dev/null +++ b/src/python_som/_accelerate.py @@ -0,0 +1,88 @@ +"""Optional numba kernel for the best-matching-unit search. Import is always safe. + +Enabled by ``pip install numba``, which this package detects and uses automatically. Without it +:func:`bmu_kernel` returns None and everything runs on the NumPy path, which stays the reference +implementation and the default. + +**Deliberately not a dependency and not an extra.** numba 0.66 requires ``numpy<2.5`` while this +package releases against 2.5, so requiring it would cap every user's NumPy; and declaring it as an +extra caps the *lockfile*, because uv resolves every extra together, which would leave development +and CI testing against an older NumPy than the release. Installing it separately keeps that +constraint in the environment that opted into it. + +**numba is imported on first use**, not when this module is imported: it costs 104 ms, and the +first training call absorbs that alongside the JIT compile. + +The kernel fuses the matrix product and the ``argmin``, keeping the running minimum in a register so +the score matrix is never written. Worth 1.0x to 2.4x, measured. Shell, not core: the kernel reaches +``_core`` as an argument rather than an import. +""" + +from __future__ import annotations + +import functools +from typing import TYPE_CHECKING + +import numpy as np + +if TYPE_CHECKING: # pragma: no cover + import numpy.typing as npt + + from ._core._protocols import BmuKernel + +__all__ = ["bmu_kernel"] + + +@functools.cache +def bmu_kernel() -> BmuKernel | None: + """Return the compiled best-matching-unit kernel, or None without the ``fast`` extra. + + Cached, so numba is imported and the kernel compiled at most once per process. + + :return: The kernel, or None. + """ + try: + from numba import njit, prange # noqa: PLC0415 deliberately deferred; see the module docs + except ImportError: + return None + + # No cover: numba compiles this, so the interpreter never executes the body and coverage cannot + # instrument it. What it does is checked by tests/test_numba_kernel.py, which asserts it returns + # the same nodes as the NumPy path. + @njit(parallel=True, cache=True) + def fused_bmu( # pragma: no cover + centred_data: npt.NDArray[np.floating], + centred_models: npt.NDArray[np.floating], + squared: npt.NDArray[np.floating], + ) -> npt.NDArray[np.intp]: + """Return the index of the nearest model for each sample, without a score matrix. + + Computes ``||w||^2 - 2 x.w`` and keeps the smallest, which orders the models exactly as + ``||x - w||`` does: the dropped ``||x||^2`` is constant per sample. Both arrays arrive + already centred, so the caller owns the cancellation fix rather than this kernel. + + ``<`` rather than ``<=``, so ties resolve to the lowest index and match ``argmin``. + + :param centred_data: Samples, shifted, of shape ``(n_samples, n_features)``. + :param centred_models: Models, shifted, of shape ``(n_nodes, n_features)``. + :param squared: Squared norm of each centred model. + :return: One flat node index per sample. + """ + n_samples, n_features = centred_data.shape + n_nodes = centred_models.shape[0] + out = np.empty(n_samples, dtype=np.intp) + for s in prange(n_samples): + best = np.inf + best_node = 0 + for node in range(n_nodes): + score = squared[node] + for f in range(n_features): + score -= 2.0 * centred_data[s, f] * centred_models[node, f] + if score < best: + best = score + best_node = node + out[s] = best_node + return out + + kernel: BmuKernel = fused_bmu # pragma: no cover only when numba is installed + return kernel # pragma: no cover diff --git a/src/python_som/_artifact.py b/src/python_som/_artifact.py index cbb1ce9..9e8253e 100644 --- a/src/python_som/_artifact.py +++ b/src/python_som/_artifact.py @@ -1,27 +1,19 @@ """Saving and loading a trained map, with the provenance needed to defend a result. Wilson et al., *Best Practices for Scientific Computing*: a result should carry its inputs, -parameters and versions. Until 0.4.0 the only way to keep a trained map was ``pickle``, which is -arbitrary code execution on load. The point here is not to forbid that but to make the safe path the -obvious one. - -**One file.** Everything lives in a single ``.npz``: the models as an array, and the metadata as a -JSON string stored alongside them. A separate sidecar was the first design and is worse, because -provenance that can be separated from its artifact will be. - -**What cannot be saved, and what happens instead.** A map holds four callables: the neighborhood, -two decays and the distance. A callable cannot be written to a file without ``pickle``, so what is -stored is its *name*, resolved on load through the registries in :mod:`python_som._core`. A map -built entirely from the shipped functions round-trips completely. One built with a caller's own -function records the name for provenance and refuses to load silently: the loader raises and names -the argument to pass it back through. - -**Security.** ``allow_pickle=False`` is passed explicitly on load, so a crafted file containing an -object array is refused by NumPy rather than executed; strategies resolve only through the -registries, so no name from the file is ever imported or evaluated; the metadata is parsed with -``json.loads``. What that buys is "cannot execute code", not "safe to load anything": an ``.npz`` is -a zip, so a hostile file can still attempt resource exhaustion through decompression. Treat one from -an untrusted source the way you would a JPEG, not the way you would a signed archive. +parameters and versions. + +**One file.** A single ``.npz`` holding the models as an array and the metadata as a JSON string. +Provenance that can be separated from its artifact will be. + +**Callables are stored by name.** A map holds four: the neighborhood, two decays and the distance. +Each is resolved on load through the registries in :mod:`python_som._core`, so a map built from the +shipped functions round-trips completely. One built with a caller's own function records the name +and refuses to load silently, naming the argument to pass it back through. + +**No pickle.** ``allow_pickle=False`` on load, so a crafted file is refused rather than executed, +and no name from a file is ever imported. See :doc:`/explanation/artifact-safety` for the limits of +that guarantee. """ from __future__ import annotations diff --git a/src/python_som/_convert.py b/src/python_som/_convert.py index bb267c7..506509f 100644 --- a/src/python_som/_convert.py +++ b/src/python_som/_convert.py @@ -1,19 +1,9 @@ """The data-input port: whatever the caller passed, in, an ``ndarray`` out. -Through 0.3.0 this module special-cased pandas, testing ``isinstance(data, pd.DataFrame | -pd.Series)`` before calling ``.to_numpy()``. That was the only use of pandas, and it was redundant: -``np.asarray`` already converts both through the ``__array__`` protocol, with identical results -including for nullable extension dtypes, which convert to ``float64`` with ``nan`` either way. - -Dropping the special case removes a required dependency and **widens** what the package -accepts, because ``__array__`` is a protocol rather than a library. polars, pyarrow, xarray and CuPy -objects all implement it and now work without python-som knowing any of them exist. Fewer -dependencies and more capability at once, which is the argument for a port rather than an adapter -per library. - -The module stays, small as it is, because it is the one place that decides what "a dataset" means. -When that decision needs to change -- a dtype policy, a shape check, an explicit error for ragged -input -- there is one place to change it, and the core keeps receiving ``ndarray`` and nothing else. +``np.asarray`` handles every input through the ``__array__`` protocol, so pandas, polars, pyarrow, +xarray and CuPy all work without this package importing any of them. + +Small, and it stays a module because it is the one place that decides what "a dataset" means. """ from __future__ import annotations diff --git a/src/python_som/_core/__init__.py b/src/python_som/_core/__init__.py index 78dd4e2..fb3f0e6 100644 --- a/src/python_som/_core/__init__.py +++ b/src/python_som/_core/__init__.py @@ -1,24 +1,18 @@ """The functional core: pure functions over NumPy arrays. -Every function here takes all of its inputs explicitly and returns a value. Nothing in this package -reads instance state, performs I/O, or knows about pandas, tqdm, or any other library beyond NumPy. -That is not a stylistic preference: it is enforced, because ruff's ``TID251`` bans those imports -everywhere except the shell modules that exist to adapt them. +Every function takes its inputs explicitly and returns a value. Nothing here reads instance state, +performs I/O, or imports anything but NumPy, which ruff's ``TID251`` enforces. -The shell around it is small by design: - -- ``python_som._convert`` converts whatever the caller passed into an ``ndarray``. The only module - that knows pandas exists. -- ``python_som._som`` holds the :class:`~python_som.SOM` class: validation, state, the training - loops, and delegation to the functions here. +The shell is ``python_som._convert``, which turns whatever the caller passed into an ``ndarray``, +``python_som._som``, which holds the state and the training loops, and ``python_som._accelerate``, +which supplies the optional kernel. References: Teuvo Kohonen, Essentials of the self-organizing map, Neural Networks 37 (2013) 52-65, https://doi.org/10.1016/j.neunet.2012.09.018 -O. J. Vrieze, Kohonen network, in: Artificial Neural Networks: An Introduction to ANN Theory and -Practice, Lecture Notes in Computer Science, vol. 931, Springer, 1995, pp. 83-100, -https://doi.org/10.1007/BFb0027024 +O. J. Vrieze, Kohonen network, in: Artificial Neural Networks, Lecture Notes in Computer Science, +vol. 931, Springer, 1995, pp. 83-100, https://doi.org/10.1007/BFb0027024 """ from __future__ import annotations @@ -31,42 +25,38 @@ ) from ._distance import euclidean_distance from ._neighborhood import ( + AXIS_PROFILES, NEIGHBORHOOD_FUNCTIONS, - NEIGHBORHOOD_KERNELS, SIGNED_NEIGHBORHOODS, + axis_matrix, axis_offsets, bubble, - bubble_kernel, + bubble_axis_profile, gaussian, - gaussian_kernel, - kernel_view, + gaussian_axis_profile, mexican_hat, - mexican_hat_kernel, - offset_span, resolve, - resolve_kernel, + resolve_axis_profile, squared_grid_distance, ) __all__ = [ + "AXIS_PROFILES", "NEIGHBORHOOD_FUNCTIONS", - "NEIGHBORHOOD_KERNELS", "SIGNED_NEIGHBORHOODS", "asymptotic_decay", + "axis_matrix", "axis_offsets", "bubble", - "bubble_kernel", + "bubble_axis_profile", "euclidean_distance", "exponential_decay", "gaussian", - "gaussian_kernel", + "gaussian_axis_profile", "inverse_decay", - "kernel_view", "linear_decay", "mexican_hat", - "mexican_hat_kernel", - "offset_span", "resolve", - "resolve_kernel", + "resolve_axis_profile", "squared_grid_distance", ] diff --git a/src/python_som/_core/_decay.py b/src/python_som/_core/_decay.py index fd82e92..b0b9eca 100644 --- a/src/python_som/_core/_decay.py +++ b/src/python_som/_core/_decay.py @@ -4,10 +4,9 @@ current value. They share one signature so that any of them, or a user-supplied equivalent, can be passed as ``learning_rate_decay`` or ``neighborhood_radius_decay``. -Kohonen (2013) does not prescribe a particular form: "The true mathematical form of sigma(t) is not -crucial, as long as its value is fairly large in the beginning of the process, say, on the order of -half of the diameter of the grid, whereafter it is gradually reduced to a fraction of it in about -1000 steps" (Section 4.1). +Kohonen (2013) Section 4.1 prescribes no particular form: "The true mathematical form of sigma(t) +is not crucial, as long as its value is fairly large in the beginning of the process ... whereafter +it is gradually reduced to a fraction of it in about 1000 steps". """ from __future__ import annotations @@ -83,12 +82,9 @@ def inverse_decay(x: float, t: int, max_t: int) -> float: } """Decay functions by name, so a saved map can name the one it used. -Each key is the function's own name, which is the least surprising mapping and the one a reader can -check against the source without a lookup table. These names are written into artifacts, so they are -public API from 0.4.0 and fixed at 1.0.0. - -A decay function is not required to be in here: a caller may pass any callable. What a name buys is -the ability to restore it from a file, and the loader says so explicitly when it cannot. +Each key is the function's own name. These are written into artifacts, so they are public API from +0.4.0 and fixed at 1.0.0. A caller may pass any callable; what a name buys is restoring it from a +file, and the loader says so when it cannot. """ diff --git a/src/python_som/_core/_linalg.py b/src/python_som/_core/_linalg.py index 5c7a4b4..f592711 100644 --- a/src/python_som/_core/_linalg.py +++ b/src/python_som/_core/_linalg.py @@ -1,17 +1,11 @@ """Principal component analysis, and the map sizing that depends on it. -Implemented on ``np.linalg.svd`` rather than scikit-learn. PCA and a z-score were the only two -things this package used scikit-learn for, and carrying it as a required dependency pulled in scipy, -joblib and threadpoolctl to reach about twenty lines of linear algebra. - -The two functions here reproduce ``sklearn.decomposition.PCA(n_components=k)`` and -``sklearn.preprocessing.StandardScaler().fit_transform`` exactly, sign conventions and -degenerate-column handling included. "Exactly" is not an assertion of faith: -``tests/test_linalg_matches_sklearn.py`` re-checks it against the real scikit-learn on every CI run, -which is why scikit-learn remains a *test* dependency. Two details are easy to get wrong and are -therefore spelled out where they are implemented: the sign convention is v-based, not the more -common u-based one, and a near-constant column is scaled by 1 rather than by its own vanishing -standard deviation. +Built on ``np.linalg.svd`` rather than scikit-learn, which it reproduces exactly, sign convention +and degenerate columns included. ``tests/test_linalg_matches_sklearn.py`` re-checks that against the +real scikit-learn on every CI run, which is why scikit-learn is still a test dependency. + +See :doc:`/explanation/why-linear-initialization-is-an-svd` for why the SVD is also the more +accurate of the two routes. """ from __future__ import annotations @@ -45,17 +39,12 @@ class PrincipalComponents(NamedTuple): def pca(data: npt.NDArray[Any], n_components: int = _N_COMPONENTS) -> PrincipalComponents: """Fit a PCA and return its mean, components and explained variance. - The singular value decomposition of the centred data gives the components directly: for - ``X - mean = U S V^T``, the rows of ``V^T`` are the component directions and ``S^2 / (n - 1)`` + For ``X - mean = U S V^T``, the rows of ``V^T`` are the component directions and ``S^2 / (n-1)`` the variance along each. - **On the sign convention.** An SVD determines each component only up to sign, so a convention is - needed for the result to be reproducible. scikit-learn's PCA calls ``svd_flip`` with - ``u_based_decision=False``, which orients each component so that its largest-magnitude *loading* - is positive. That is the less common of the two conventions in that helper, and taking its - default instead would flip the sign of some components: the fit would still be a valid PCA, but - it would not be the same one, and linear initialization would lay its models out reversed along - that axis. + The sign convention is **v-based**: each component is oriented so its largest-magnitude loading + is positive, matching scikit-learn's ``svd_flip(..., u_based_decision=False)``. The other + convention would give a valid but different PCA, and lay the initial models out reversed. :param data: Array of shape ``(n_samples, n_features)``. :param n_components: Number of components to keep. diff --git a/src/python_som/_core/_maps.py b/src/python_som/_core/_maps.py index 2b33796..25816aa 100644 --- a/src/python_som/_core/_maps.py +++ b/src/python_som/_core/_maps.py @@ -7,7 +7,7 @@ import numpy as np -from ._match import winner +from ._match import bmu_indices, winner from ._neighborhood import bubble if TYPE_CHECKING: # pragma: no cover @@ -31,17 +31,12 @@ def u_matrix( ) -> npt.NDArray[np.floating]: """Return the U-matrix: the summed distance from each model to its immediate neighbours. - Ultsch's display (1993), cited by Kohonen (2013) Section 3.6 as the way cluster structure is - made visible on the grid: a large value means neighbouring models are far apart, so it reads - as a boundary. + Ultsch's display (1993), cited by Kohonen (2013) Section 3.6: a large value means neighbouring + models are far apart, so it reads as a boundary. - The adjacency is deliberately a flat ring of radius 1 rather than the configured neighborhood - function. The U-matrix describes the grid, not the training schedule. The centre is included and - contributes a distance of zero, so it does not affect the sum. - - Distances are computed and consumed one node at a time rather than accumulated into a full - ``(x, y, x, y)`` tensor, which would cost ``(x*y)**2`` floats: about 800 MB on a 100x100 map, to - produce ``x*y`` numbers. + The adjacency is a flat ring of radius 1 rather than the configured neighborhood, because the + U-matrix describes the grid and not the training schedule. Distances are consumed one node at a + time; the full ``(x, y, x, y)`` tensor would be 800 MB on a 100x100 map. :param weights: Models, of shape ``(x, y, n_features)``. :param shape: Shape of the grid. @@ -105,8 +100,9 @@ def winner_map( result: dict[tuple[int, int], list[npt.NDArray[Any]]] = { (int(i), int(j)): [] for i, j in np.ndindex(shape) } - for sample in data: - result[winner(sample, weights, distance)].append(sample) + rows, columns = np.unravel_index(bmu_indices(data, weights, distance), shape) + for sample, row, column in zip(data, rows, columns, strict=True): + result[int(row), int(column)].append(sample) return result @@ -136,6 +132,7 @@ def label_map( counts: dict[tuple[int, int], Counter[Any]] = { (int(i), int(j)): Counter() for i, j in np.ndindex(shape) } - for sample, label in zip(data, labels, strict=True): - counts[winner(sample, weights, distance)].update([label]) + rows, columns = np.unravel_index(bmu_indices(data, weights, distance), shape) + for label, row, column in zip(labels, rows, columns, strict=True): + counts[int(row), int(column)].update([label]) return counts diff --git a/src/python_som/_core/_match.py b/src/python_som/_core/_match.py index fcaf46a..890a95d 100644 --- a/src/python_som/_core/_match.py +++ b/src/python_som/_core/_match.py @@ -10,12 +10,19 @@ import numpy as np +from ._distance import euclidean_distance + if TYPE_CHECKING: # pragma: no cover import numpy.typing as npt - from ._protocols import DistanceFunction + from ._protocols import BmuKernel, DistanceFunction + +__all__ = ["accumulate", "activate", "bmu_indices", "quantization", "winner"] -__all__ = ["accumulate", "activate", "quantization", "winner"] +#: Bytes the winner search may hold in its score block at once, setting the chunk size. Tuned: at +#: 60x60 with 2000 samples an 8 MB budget is 2.6x slower and 8x heavier, because a block that fits +#: in cache is read back by ``argmin`` for free. See /explanation/how-batch-training-is-computed. +_SCORE_BUDGET_BYTES = 512_000 def activate( @@ -59,7 +66,67 @@ def quantization( :param distance: Dissimilarity measure. :return: One distance per sample. """ - return np.array([distance(i, weights[winner(i, weights, distance)]) for i in data]) + flat = weights.reshape(-1, weights.shape[-1]) + nodes = bmu_indices(data, weights, distance) + # The distance is recomputed against the chosen model rather than read out of the search, which + # keeps this exact for the Euclidean case: `bmu_indices` drops ||x||^2, so its scores order the + # models correctly but are not distances. + return np.array([distance(x, flat[node]) for x, node in zip(data, nodes, strict=True)]) + + +def bmu_indices( + data: npt.NDArray[Any], + weights: npt.NDArray[Any], + distance: DistanceFunction, + kernel: BmuKernel | None = None, +) -> npt.NDArray[np.intp]: + """Return the flat index of the best-matching model for every sample. + + This is Eq. (4) of Kohonen (2013), ``c = argmin_i ||x - m_i||``, for a whole dataset. Ties go to + the first index in C order, which is ``argmin``'s behaviour and matches :func:`winner`. + + For the Euclidean distance this expands the norm and drops the ``||x||^2`` term, which is + constant across models, leaving a matrix product. Any other distance takes the loop, since only + the Euclidean one has that identity. + + **Not the dot-product map of Kohonen Section 4.5**, which is a different algorithm requiring + renormalized models. This is an exact re-expansion of the Euclidean distance. + + **The centring is not an optimization.** Without it the expansion cancels catastrophically: + with models offset by 1e9, 499 of 500 samples get a different node. See + :doc:`/explanation/how-batch-training-is-computed`. + + :param data: Dataset of shape ``(n_samples, n_features)``. + :param weights: Models, of shape ``(x, y, n_features)``. + :param distance: Dissimilarity measure. + :param kernel: Optional accelerated search, from ``python_som._accelerate``. Passed in rather + than imported, so this module stays numpy-only. + :return: One flat node index per sample. + """ + flat = weights.reshape(-1, weights.shape[-1]) + if distance is not euclidean_distance: + return np.array([np.asarray(distance(x, flat)).argmin() for x in data], dtype=np.intp) + + shift = flat.mean(axis=0) + centred = flat - shift + squared = np.einsum("nf,nf->n", centred, centred) + + if kernel is not None: # pragma: no cover reached only when numba is installed + return kernel(data - shift, centred, squared) + + n_nodes = len(flat) + chunk = max(1, _SCORE_BUDGET_BYTES // (n_nodes * 8)) + scores = np.empty((chunk, n_nodes)) + out = np.empty(len(data), dtype=np.intp) + for start in range(0, len(data), chunk): + block = data[start : start + chunk] + # Into a preallocated buffer: allocating one per chunk was 1.5x slower and 8x heavier. + np.matmul(block - shift, centred.T, out=scores[: len(block)]) + block_scores = scores[: len(block)] + block_scores *= -2.0 + block_scores += squared + out[start : start + len(block)] = block_scores.argmin(axis=1) + return out def accumulate( @@ -67,6 +134,7 @@ def accumulate( weights: npt.NDArray[Any], shape: tuple[int, int], distance: DistanceFunction, + kernel: BmuKernel | None = None, ) -> tuple[npt.NDArray[np.floating], npt.NDArray[np.floating]]: """Sum the samples mapped to each node, and count them. @@ -77,12 +145,12 @@ def accumulate( :param weights: Models, of shape ``(x, y, n_features)``. :param shape: Shape of the grid. :param distance: Dissimilarity measure. + :param kernel: Optional accelerated search; see :func:`bmu_indices`. :return: Per-node sums of shape ``(x, y, n_features)`` and counts of shape ``(x, y)``. """ - sums = np.zeros((*shape, weights.shape[-1])) - counts = np.zeros(shape) - for sample in data: - node = winner(sample, weights, distance) - sums[node] += sample - counts[node] += 1 - return sums, counts + nodes = bmu_indices(data, weights, distance, kernel) + n_nodes = shape[0] * shape[1] + sums = np.zeros((n_nodes, weights.shape[-1])) + np.add.at(sums, nodes, data) + counts = np.bincount(nodes, minlength=n_nodes).astype(float) + return sums.reshape(*shape, weights.shape[-1]), counts.reshape(shape) diff --git a/src/python_som/_core/_neighborhood.py b/src/python_som/_core/_neighborhood.py index 662f9d1..8154b65 100644 --- a/src/python_som/_core/_neighborhood.py +++ b/src/python_som/_core/_neighborhood.py @@ -1,24 +1,16 @@ """Neighborhood functions: how the winner's correction spreads over the grid. -Kohonen (2013) Eq. (5) defines the neighborhood as a function of ``sqdist(c, i)``, "the square of -the geometric distance between the nodes c and i in the grid". Vrieze (1995) Fig. 3 plots the -"Mexican-hat" lateral interaction against a single axis labelled "Lateral distance", writes the -coefficient as ``h_{i i_c} = 1 / ||i_c - i||``, and states that the grid is assumed to be a metric -space. - -The consequence is that a neighborhood function must depend on the distance between two nodes and -not on the two axis offsets separately. The gaussian happens to factor into a product of per-axis -terms, but that is a property of the exponential, not of neighborhood functions in general: an -outer product of two 1-D Ricker wavelets is positive in the diagonal quadrants where both factors -are negative, placing an excitatory lobe exactly where the mexican hat must inhibit. +Kohonen (2013) Eq. (5) defines a neighborhood as a function of ``sqdist(c, i)``, the squared grid +distance between two nodes, so it must depend on that distance and not on the two axis offsets +separately. See :doc:`/explanation/why-isotropy-matters` for why an outer product of two 1-D +profiles is wrong for anything but the gaussian. References: Teuvo Kohonen, Essentials of the self-organizing map, Neural Networks 37 (2013) 52-65, https://doi.org/10.1016/j.neunet.2012.09.018 -O. J. Vrieze, Kohonen network, in: Artificial Neural Networks: An Introduction to ANN Theory and -Practice, Lecture Notes in Computer Science, vol. 931, Springer, 1995, pp. 83-100, -https://doi.org/10.1007/BFb0027024 +O. J. Vrieze, Kohonen network, in: Artificial Neural Networks, Lecture Notes in Computer Science, +vol. 931, Springer, 1995, pp. 83-100, https://doi.org/10.1007/BFb0027024 """ from __future__ import annotations @@ -28,25 +20,24 @@ import numpy as np import numpy.typing as npt -from ._protocols import KernelFunction, NeighborhoodFunction +from ._protocols import AxisProfile, KernelFunction, NeighborhoodFunction __all__ = [ + "AXIS_PROFILES", "NEIGHBORHOOD_FUNCTIONS", - "NEIGHBORHOOD_KERNELS", "SIGNED_NEIGHBORHOODS", + "AxisProfile", "KernelFunction", "NeighborhoodFunction", + "axis_matrix", "axis_offsets", "bubble", - "bubble_kernel", + "bubble_axis_profile", "gaussian", - "gaussian_kernel", - "kernel_view", + "gaussian_axis_profile", "mexican_hat", - "mexican_hat_kernel", - "offset_span", "resolve", - "resolve_kernel", + "resolve_axis_profile", "squared_grid_distance", ] @@ -73,9 +64,8 @@ def _validate_radius(sigma: float, *, allow_zero: bool = False) -> None: def axis_offsets(length: int, center: int, *, cyclic: bool) -> npt.NDArray[np.floating]: """Signed offsets from ``center`` to every coordinate along one axis. - On a cyclic axis the minimum-image convention folds each offset into ``[-length/2, length/2)``, - so the shortest way around the torus is used. Both tails must be folded: an offset of -9 on an - axis of length 10 represents a distance of 1, not 9. + On a cyclic axis the minimum-image convention folds each offset into ``[-length/2, length/2)``. + Both tails must be folded: an offset of -9 on an axis of length 10 is a distance of 1, not 9. :param length: Number of nodes along the axis. :param center: Coordinate of the winner along the axis. @@ -88,28 +78,6 @@ def axis_offsets(length: int, center: int, *, cyclic: bool) -> npt.NDArray[np.fl return d -def offset_span(length: int, *, cyclic: bool) -> npt.NDArray[np.floating]: - """Every offset any pair of nodes on this axis can have: ``-(length-1) .. (length-1)``. - - The full-range counterpart of :func:`axis_offsets`, which gives the offsets from one particular - centre. Because a neighborhood depends on the offset alone and never on where the winner sits, - one array over this span serves every node -- which is what lets batch training evaluate the - neighborhood once per iteration instead of once per node. - - The cyclic fold is the same minimum-image convention, applied with the real period ``length`` - rather than the span's own width. That is the whole reason this cannot be expressed as - ``axis_offsets`` on a ``2*length-1`` axis: the fold would then use the wrong period. - - :param length: Number of nodes along the axis. - :param cyclic: Whether the axis wraps around. - :return: Signed offsets, ``2 * length - 1`` of them, centred on zero. - """ - d = np.arange(-(length - 1), length, dtype=float) - if cyclic: - d = (d + length / 2) % length - length / 2 - return d - - def squared_grid_distance( shape: Grid, c: Coordinates, cyclic: tuple[bool, bool] ) -> npt.NDArray[np.floating]: @@ -125,18 +93,8 @@ def squared_grid_distance( return np.add.outer(np.square(dx), np.square(dy)) -# --------------------------------------------------------------------------------------------- -# The profiles: one implementation of each formula, shared by the per-node and kernel forms. -# -# Each takes the two axes' offsets rather than a grid and a centre, because that is the only thing -# the two forms differ in: the per-node function passes `axis_offsets` from one winner, and the -# kernel builder passes `offset_span` covering every winner at once. Keeping one copy of the formula -# is what makes the two bit-identical by construction rather than by agreement, which matters here: -# the defect that started this whole investigation was a plausible-looking second version of the -# mexican hat that disagreed with the first. -# -# Validation lives here, so neither form can skip it. -# --------------------------------------------------------------------------------------------- +# One implementation of each formula, taking axis offsets rather than a grid and a centre. +# Validation lives here, so no caller can skip it. def _gaussian_profile( @@ -175,10 +133,8 @@ def _bubble_profile( ) -> npt.NDArray[np.floating]: """Evaluate the Chebyshev indicator ``max(|dx|, |dy|) <= round(sigma)``. - Deliberately not built on ``sqdist``, unlike the other two: the bubble's metric is Chebyshev, so - it is a product of two per-axis indicators rather than a function of a Euclidean distance. That - asymmetry is the implementation following Vrieze's appendix, and is preserved rather than - quietly unified. See :func:`bubble`. + Not built on ``sqdist``, unlike the other two: the bubble's metric is Chebyshev. See + :func:`bubble`. :param dx: Offsets along the first axis. :param dy: Offsets along the second axis. @@ -196,8 +152,8 @@ def gaussian( ) -> npt.NDArray[np.floating]: """Gaussian neighborhood, ``exp(-sqdist(c, i) / (2 * sigma**2))``. - This is Eq. (5) of Kohonen (2013) with the learning rate factored out, so ``h(c, c) == 1``. - Strictly positive everywhere and monotonically decreasing with distance. + Eq. (5) of Kohonen (2013) with the learning rate factored out, so ``h(c, c) == 1``. Strictly + positive and monotonically decreasing with distance. :param shape: Shape of the network. :param c: Coordinates of the winner. @@ -218,19 +174,12 @@ def mexican_hat( ) -> npt.NDArray[np.floating]: """Mexican hat neighborhood, ``(1 - u) * exp(-u)`` over ``u = sqdist(c, i) / (2 * sigma**2)``. - Also known as the Ricker wavelet or the Laplacian of Gaussian. This is the biologically - motivated lateral-interaction function: nodes near the winner are excited, nodes past a certain - distance are inhibited, and the inhibition vanishes as distance grows further (Vrieze 1995, - Fig. 3). + The Ricker wavelet, or Laplacian of Gaussian: excitatory near the winner, inhibitory beyond it, + vanishing with distance (Vrieze 1995, Fig. 3). Normalized so ``h(c, c) == 1``, zero at + ``sqrt(2) * sigma``, minimum ``-exp(-2)`` at ``2 * sigma``. - Normalized so ``h(c, c) == 1``. Crosses zero at a radius of ``sqrt(2) * sigma`` and reaches its - minimum of ``-exp(-2)``, about -0.135, at a radius of ``2 * sigma``. - - This is deliberately not the outer product of two 1-D Ricker wavelets. See the module docstring - for why that construction is wrong. - - Takes negative values, so it cannot be used with batch training; see - :data:`SIGNED_NEIGHBORHOODS`. + Not an outer product of two 1-D Ricker wavelets, which is a different and wrong function; see + :doc:`/explanation/why-isotropy-matters`. Signed, so batch training rejects it. :param shape: Shape of the network. :param c: Coordinates of the winner. @@ -251,27 +200,17 @@ def bubble( ) -> npt.NDArray[np.floating]: """Flat neighborhood: 1 for nodes within ``sigma`` of the winner, 0 elsewhere. - This is the truncated inner, excitatory lobe of the mexican hat. Vrieze (1995) p. 85: "In - Kohonen networks usually only the inner stimulation area is used, i.e., when a neuron i fires, a - positive feedback takes place for all neurons i', whose distance to i is smaller than some given - number rho", and notes that this flat choice is "just as effective and sometimes even better" - than a distance-dependent one. - - **The metric here is Chebyshev, not Euclidean**, so the region is a square rather than a disc: - a node is included when ``max(|dx|, |dy|) <= round(sigma)``. That matches the pseudo-code in - Vrieze's appendix, which computes ``b = MAX(ABS(i - w_i), ABS(j - w_j))``, though Kohonen's - phrase "up to a certain radius from the winner" (Section 4.1) reads as Euclidean. The two - sources genuinely differ; this implementation follows Vrieze, and the choice is preserved rather - than changed so that existing results stay reproducible. + The truncated inner lobe of the mexican hat, which Vrieze (1995) p. 85 calls "just as effective + and sometimes even better" than a distance-dependent one. - One consequence worth stating, because it is easy to assume otherwise: a Chebyshev ball is not - isotropic under the Euclidean metric. On a large enough grid, nodes at equal Euclidean distance - from the winner can fall on opposite sides of the boundary. The smallest case is a radius of - ``sqrt(50)``, where ``(5, 5)`` lies inside a ``sigma = 5`` square while ``(7, 1)`` lies outside. + **The metric is Chebyshev, not Euclidean**, so the region is a square: a node is included when + ``max(|dx|, |dy|) <= round(sigma)``. This follows Vrieze's appendix, where Kohonen's "up to a + certain radius" (Section 4.1) reads as Euclidean; the two sources differ, and + :doc:`/explanation/why-isotropy-matters` covers the consequence, that a Chebyshev ball is not + isotropic under the Euclidean metric. - Unlike the other neighborhood functions, a radius of zero is admissible: it selects the winner - alone, which is well defined, whereas for the gaussian and the mexican hat it is a division by - zero. + A radius of zero is admissible here and not for the other two: it selects the winner alone, + where they would divide by zero. :param shape: Shape of the network. :param c: Coordinates of the winner. @@ -296,126 +235,107 @@ def bubble( } """Neighborhood functions by name. ``mexican_hat`` is an alias of ``mexicanhat``.""" -SIGNED_NEIGHBORHOODS: Final[frozenset[str]] = frozenset({"mexicanhat", "mexican_hat"}) -"""Names of neighborhood functions that take negative values, which batch training cannot use.""" - - -def resolve(name: str) -> NeighborhoodFunction: - """Look up a neighborhood function by name. - - :param name: Name of the neighborhood function. - :return: The corresponding function. - :raises ValueError: If the name is not recognised. - """ - try: - return NEIGHBORHOOD_FUNCTIONS[name] - except KeyError as exc: - valid = sorted(NEIGHBORHOOD_FUNCTIONS) - msg = ( - f"Invalid value for 'neighborhood_function' parameter: {name!r}. " - f"Value should be one of {valid}" - ) - raise ValueError(msg) from exc +# Axis profiles: the per-axis factor of a separable neighborhood, used to contract Eq. (8) into two +# matrix products. A contraction strategy for a function of sqdist, never a redefinition of one. +# The mexican hat has no entry and must not acquire one: (1 - u) exp(-u) does not factor. +# See /explanation/how-batch-training-is-computed. -# --------------------------------------------------------------------------------------------- -# Kernels: one evaluation per iteration instead of one per node. -# --------------------------------------------------------------------------------------------- +def gaussian_axis_profile(d: npt.NDArray[np.floating], sigma: float) -> npt.NDArray[np.floating]: + """Per-axis factor of the gaussian, ``exp(-d^2 / (2 sigma^2))``. -def gaussian_kernel( - shape: Grid, sigma: float, cyclic: tuple[bool, bool] -) -> npt.NDArray[np.floating]: - """Evaluate the gaussian over every offset, to be sliced per node by :func:`kernel_view`. + Its product over the two axes is :func:`gaussian`, because the exponential factors. - :param shape: Shape of the network. - :param sigma: Neighborhood radius. - :param cyclic: Whether each axis wraps around. - :return: Weights of shape ``(2 * shape[0] - 1, 2 * shape[1] - 1)``. + :param d: Offsets along one axis. + :param sigma: Neighborhood radius. Must be finite and positive. + :return: Weights for those offsets. :raises ValueError: If the radius is not a finite positive number. """ - return _gaussian_profile( - offset_span(shape[0], cyclic=cyclic[0]), offset_span(shape[1], cyclic=cyclic[1]), sigma - ) - - -def mexican_hat_kernel( - shape: Grid, sigma: float, cyclic: tuple[bool, bool] -) -> npt.NDArray[np.floating]: - """Evaluate the mexican hat over every offset. See :func:`gaussian_kernel`. + _validate_radius(sigma) + return np.exp(-np.square(d) / (2.0 * sigma * sigma)) - :param shape: Shape of the network. - :param sigma: Neighborhood radius. - :param cyclic: Whether each axis wraps around. - :return: Weights of shape ``(2 * shape[0] - 1, 2 * shape[1] - 1)``. - :raises ValueError: If the radius is not a finite positive number. - """ - return _mexican_hat_profile( - offset_span(shape[0], cyclic=cyclic[0]), offset_span(shape[1], cyclic=cyclic[1]), sigma - ) +def bubble_axis_profile(d: npt.NDArray[np.floating], sigma: float) -> npt.NDArray[np.floating]: + """Per-axis factor of the bubble, the indicator ``|d| <= round(sigma)``. -def bubble_kernel(shape: Grid, sigma: float, cyclic: tuple[bool, bool]) -> npt.NDArray[np.floating]: - """Evaluate the bubble over every offset. See :func:`gaussian_kernel`. + Its product over the two axes is :func:`bubble`. It factors because the metric is Chebyshev; a + Euclidean disc would not. - :param shape: Shape of the network. - :param sigma: Neighborhood radius. - :param cyclic: Whether each axis wraps around. - :return: Weights of shape ``(2 * shape[0] - 1, 2 * shape[1] - 1)``. + :param d: Offsets along one axis. + :param sigma: Neighborhood radius, rounded to the nearest integer. Must be finite and + non-negative. + :return: Weights for those offsets. :raises ValueError: If the radius is not a finite non-negative number. """ - return _bubble_profile( - offset_span(shape[0], cyclic=cyclic[0]), offset_span(shape[1], cyclic=cyclic[1]), sigma - ) + _validate_radius(sigma, allow_zero=True) + return (np.abs(d) <= int(np.around(sigma))).astype(float) -NEIGHBORHOOD_KERNELS: Final[dict[str, KernelFunction]] = { - "gaussian": gaussian_kernel, - "bubble": bubble_kernel, - "mexicanhat": mexican_hat_kernel, - "mexican_hat": mexican_hat_kernel, +AXIS_PROFILES: Final[dict[str, AxisProfile]] = { + "gaussian": gaussian_axis_profile, + "bubble": bubble_axis_profile, } -"""Kernel form of each neighborhood function, keyed exactly as :data:`NEIGHBORHOOD_FUNCTIONS`. +"""Per-axis factor of each separable neighborhood, keyed as :data:`NEIGHBORHOOD_FUNCTIONS`. -Every registered name has one, which is what lets batch training use the kernel path unconditionally -instead of carrying a fallback branch for a case that cannot arise. +Membership decides what batch training accepts: a neighborhood absent from it is rejected by name +rather than approximated. """ -def kernel_view( - kernel: npt.NDArray[np.floating], shape: Grid, c: Coordinates -) -> npt.NDArray[np.floating]: - """Extract node ``c``'s neighborhood from a kernel, as a view rather than a copy. +def resolve_axis_profile(name: str) -> AxisProfile: + """Look up the per-axis factor of a neighborhood function by name. - The kernel is indexed by offset, with offset zero at ``(shape[0] - 1, shape[1] - 1)``. Node - ``c`` sees offsets ``i - c`` for each node ``i``, so its neighborhood is the ``shape``-sized - block starting at ``(shape[0] - 1 - c[0], shape[1] - 1 - c[1])``. + :param name: Name of the neighborhood function. + :return: The corresponding axis profile. + :raises ValueError: If the function has no axis profile, and so is not separable. + """ + try: + return AXIS_PROFILES[name] + except KeyError as exc: + valid = sorted(AXIS_PROFILES) + msg = ( + f"The {name!r} neighborhood function is not separable, so it has no axis profile. " + f"Value should be one of {valid}" + ) + raise ValueError(msg) from exc - Returning a view is the point: copying ``shape[0] * shape[1]`` floats per node would give back - much of what evaluating the kernel once saved. Downstream reads it only, and ``np.sum`` and - ``np.einsum`` are both happy with a non-contiguous view. - :param kernel: Kernel from one of the ``*_kernel`` functions. - :param shape: Shape of the network. - :param c: Coordinates of the node whose neighborhood is wanted. - :return: A read-only-by-convention view of shape ``shape``. +def axis_matrix( + length: int, sigma: float, *, cyclic: bool, profile: AxisProfile +) -> npt.NDArray[np.floating]: + """Build ``H[a, c] = profile(a - c)`` for every pair of coordinates on one axis. + + Contracting ``sums`` against one of these per axis evaluates Eq. (8) for every node at once. The + cyclic fold is the minimum-image convention of :func:`axis_offsets`, on pairwise offsets. + + :param length: Number of nodes along the axis. + :param sigma: Neighborhood radius. + :param cyclic: Whether the axis wraps around. + :param profile: Per-axis factor to evaluate, from :data:`AXIS_PROFILES`. + :return: Weights of shape ``(length, length)``. """ - return kernel[ - shape[0] - 1 - c[0] : 2 * shape[0] - 1 - c[0], shape[1] - 1 - c[1] : 2 * shape[1] - 1 - c[1] - ] + d = np.subtract.outer(np.arange(length), np.arange(length)).astype(float) + if cyclic: + d = (d + length / 2) % length - length / 2 + return profile(d, sigma) + + +SIGNED_NEIGHBORHOODS: Final[frozenset[str]] = frozenset({"mexicanhat", "mexican_hat"}) +"""Names of neighborhood functions that take negative values, which batch training cannot use.""" -def resolve_kernel(name: str) -> KernelFunction: - """Look up the kernel form of a neighborhood function by name. +def resolve(name: str) -> NeighborhoodFunction: + """Look up a neighborhood function by name. :param name: Name of the neighborhood function. - :return: The corresponding kernel builder. + :return: The corresponding function. :raises ValueError: If the name is not recognised. """ try: - return NEIGHBORHOOD_KERNELS[name] + return NEIGHBORHOOD_FUNCTIONS[name] except KeyError as exc: - valid = sorted(NEIGHBORHOOD_KERNELS) + valid = sorted(NEIGHBORHOOD_FUNCTIONS) msg = ( f"Invalid value for 'neighborhood_function' parameter: {name!r}. " f"Value should be one of {valid}" diff --git a/src/python_som/_core/_protocols.py b/src/python_som/_core/_protocols.py index 785facc..01742d4 100644 --- a/src/python_som/_core/_protocols.py +++ b/src/python_som/_core/_protocols.py @@ -1,18 +1,14 @@ -"""Contracts for the three strategies a caller can replace. +"""Contracts for the strategies a caller can replace. -A neighborhood, a decay and a distance are all things a user may supply their own version of. Typed -as bare ``Callable[...]`` aliases, mypy checks little more than the argument count; as Protocols it -checks the shape of the call against a named contract, and the error names the protocol rather than -printing two structural types side by side. +Protocols rather than bare ``Callable`` aliases, so mypy checks the shape of the call against a +named contract instead of only the argument count. -**Every parameter is positional-only** (the ``/`` in each ``__call__``). Without it, a Protocol -requires the *names* to match as well as the types, so a user's ``def my_decay(rate, step, total)`` -would fail against a protocol that named them differently. Positional-only says what is actually -true: these are called positionally, and only the order and the types matter. +**Every parameter is positional-only.** Without the ``/``, a Protocol also requires the parameter +*names* to match, so a user's ``def my_decay(rate, step, total)`` would fail against a protocol that +named them differently. -These are structural, so nothing needs to inherit from them. Every function already in the package -satisfies its protocol, and so does any existing user-supplied callable with the right signature -- -this adds checking, not a requirement. +Structural, so nothing inherits from them: this adds checking, not a requirement. See +:doc:`/how-to/use-a-custom-strategy`. """ from __future__ import annotations @@ -23,16 +19,23 @@ import numpy as np import numpy.typing as npt -__all__ = ["DecayFunction", "DistanceFunction", "KernelFunction", "NeighborhoodFunction"] +__all__ = [ + "AxisProfile", + "BmuKernel", + "DecayFunction", + "DistanceFunction", + "KernelFunction", + "NeighborhoodFunction", +] @runtime_checkable class NeighborhoodFunction(Protocol): """Weights the winner's correction across the grid, as a function of grid distance. - Kohonen (2013) Eq. (5) requires this to depend on ``sqdist(c, i)`` alone -- the distance between - two nodes -- not on the two axis offsets separately. A separable product of per-axis profiles - satisfies the signature but is only correct for the gaussian. + Kohonen (2013) Eq. (5) requires this to depend on ``sqdist(c, i)`` alone, not on the two axis + offsets separately. A separable product satisfies the signature and is only correct for the + gaussian. """ def __call__( @@ -74,7 +77,7 @@ class DistanceFunction(Protocol): """Dissimilarity between an input vector and one or many models. Called both with a single model and with the whole ``(x, y, n_features)`` array, so an - implementation must broadcast over leading axes rather than assume one vector. + implementation must broadcast over leading axes. """ def __call__(self, x: Any, weights: Any, /) -> npt.NDArray[np.floating]: # noqa: ANN401 @@ -87,12 +90,61 @@ def __call__(self, x: Any, weights: Any, /) -> npt.NDArray[np.floating]: # noqa ... +@runtime_checkable +class BmuKernel(Protocol): + """An accelerated best-matching-unit search, supplied from outside the core. + + Optional. ``python_som._accelerate`` provides one with the ``fast`` extra installed; otherwise + the NumPy path in :func:`~python_som._core._match.bmu_indices` runs. Passed as an argument + rather than imported, so the core stays numpy-only. + + Both arrays arrive already shifted by a common vector, which is what stops the expanded norm + cancelling far from the origin. A kernel must not shift them again. + """ + + def __call__( + self, + centred_data: npt.NDArray[np.floating], + centred_models: npt.NDArray[np.floating], + squared: npt.NDArray[np.floating], + /, + ) -> npt.NDArray[np.intp]: + """Return the index of the nearest model for each sample. + + :param centred_data: Samples, shifted, of shape ``(n_samples, n_features)``. + :param centred_models: Models, shifted, of shape ``(n_nodes, n_features)``. + :param squared: Squared norm of each centred model. + :return: One flat node index per sample, ties going to the lowest index. + """ + ... + + +@runtime_checkable +class AxisProfile(Protocol): + """The per-axis factor of a separable neighborhood, over offsets along one axis. + + Defined only where the factorisation is an identity, which is the gaussian and the bubble. Not a + general way to build a neighborhood: :class:`NeighborhoodFunction` remains the definition. + """ + + def __call__(self, d: npt.NDArray[np.floating], sigma: float, /) -> npt.NDArray[np.floating]: + """Evaluate the factor over offsets along one axis. + + :param d: Offsets along the axis. + :param sigma: Neighborhood radius. + :return: Weights for those offsets. + """ + ... + + @runtime_checkable class KernelFunction(Protocol): """A neighborhood evaluated over every offset at once, independent of any particular winner. - The kernel form of a :class:`NeighborhoodFunction`, used by batch training so that the - neighborhood is computed once per iteration rather than once per node. + .. deprecated:: 0.7.0 + Batch training now contracts an :class:`AxisProfile` per axis, and nothing in the package + produces a kernel. Retained because it is part of the public surface; it will be removed at + 1.0.0. """ def __call__( diff --git a/src/python_som/_core/_update.py b/src/python_som/_core/_update.py index 765baaf..421182b 100644 --- a/src/python_som/_core/_update.py +++ b/src/python_som/_core/_update.py @@ -1,24 +1,8 @@ """The two update rules, as pure functions returning new models. -Both return a new array rather than mutating their argument, which is what lets them be tested -without constructing a :class:`~python_som.SOM`, and what 1.0.0's ``fit`` will assign to -``weights_``. - -That choice costs a little speed rather than gaining it, which is worth stating plainly because an -earlier draft of this module claimed the opposite. Measured against the in-place -``weights += alpha * h[..., None] * (sample - weights)`` that 0.3.0 shipped, with the two arms -interleaved and compared on medians: the pure form is **roughly 10% slower on small maps** (20x20, -50x50) and **indistinguishable on large ones** (100x100 and up, where the interquartile ranges -overlap). It is never faster. On a 20x20 map the penalty is single-digit milliseconds across a -10,000-iteration run, which is not a reason to give up a function that can be tested without -constructing a network. - -The claim it replaces was that the pure form ran up to 2.9x *faster*. That came from a benchmark -whose two arms did not compute quite the same thing and whose repeats were not interleaved, so -thermal drift was read as a speedup; it does not replicate. ``benchmarks/bench_update.py`` is the -corrected version, and it asserts the two forms agree at exactly ``0.0`` before it will report a -timing. Run it rather than trusting the summary above, since the ratios depend on the machine. The -equality itself is a test, in ``tests/test_core_boundary.py``. +Both return a new array rather than mutating their argument, so they can be tested without +constructing a :class:`~python_som.SOM`. That costs roughly 10% on small maps and nothing on large +ones; ``benchmarks/bench_update.py`` measures it. """ from __future__ import annotations @@ -28,13 +12,8 @@ import numpy as np if TYPE_CHECKING: # pragma: no cover - from collections.abc import Callable - import numpy.typing as npt - #: Given node coordinates, return that node's neighborhood over the grid. - NeighborhoodOf = Callable[[tuple[int, int]], npt.NDArray[np.floating]] - __all__ = ["batch_update", "stepwise_update"] @@ -69,39 +48,36 @@ def batch_update( weights: npt.NDArray[Any], sums: npt.NDArray[np.floating], counts: npt.NDArray[np.floating], - neighborhood_of: NeighborhoodOf, - shape: tuple[int, int], + hx: npt.NDArray[np.floating], + hy: npt.NDArray[np.floating], ) -> npt.NDArray[np.floating]: """Recompute every model as the neighborhood-weighted mean of the data around it. This is Eq. (8) of Kohonen (2013), ``m_i = sum_j n_j h_ji xbar_j / sum_j n_j h_ji``, where ``sums[j]`` is ``n_j * xbar_j``. - Two properties are worth stating because they are easy to get wrong: + The sum over node pairs is a convolution, and a separable ``h`` contracts it to two matrix + products with no loop over nodes. See :doc:`/explanation/how-batch-training-is-computed`. - **A model with no data in its neighborhood keeps its previous value.** Building the result from - a zeroed array instead destroys it; on a 30x30 map with 20 samples and a small radius that - wiped 282 of 900 models in a single step. + Three invariants: - **The denominator needs no tolerance, only ``> 0``.** Every term of ``sum_j n_j h_ji`` is - non-negative, because a signed neighborhood cannot reach this function: batch training rejects - the mexican hat, and a caller cannot supply an arbitrary neighborhood since only registered - names resolve. A sum of non-negative floats admits no cancellation, so it is zero exactly when - every term is zero, which is exactly the "no data in reach" case. An epsilon here would be an - invented number guarding a condition that cannot arise. + - Every model is computed from the models as they stood at the start of the iteration, which is + the concurrent update Kohonen requires in Section 4.4. + - A model with no data in its neighborhood keeps its previous value. Building the result from a + zeroed array wiped 282 of 900 models in one step on a 30x30 map; ``out=`` with ``where=`` is + what preserves it. + - The denominator needs no tolerance, only ``> 0``. Every term is non-negative, since batch + training rejects signed neighborhoods, so the sum is zero exactly when every term is. :param weights: Current models, of shape ``(x, y, n_features)``. :param sums: Per-node sums of the samples mapped to each node. :param counts: Per-node counts of the samples mapped to each node. - :param neighborhood_of: Callable taking node coordinates and returning its neighborhood. - :param shape: Shape of the grid. + :param hx: Per-axis neighborhood factor for the first axis, of shape ``(x, x)``. + :param hy: Per-axis neighborhood factor for the second axis, of shape ``(y, y)``. :return: The updated models, as a new array. """ + numerator = np.einsum("ac,bd,cdf->abf", hx, hy, sums, optimize=True) + denominator = np.einsum("ac,bd,cd->ab", hx, hy, counts, optimize=True) updated = weights.copy() - for node in np.ndindex(shape): - node_2d = (int(node[0]), int(node[1])) - h = neighborhood_of(node_2d) - denominator = float(np.sum(h * counts)) - if denominator > 0: - updated[node_2d] = np.einsum("xy,xyf->f", h, sums) / denominator + np.divide(numerator, denominator[..., None], out=updated, where=denominator[..., None] > 0) return updated diff --git a/src/python_som/_enums.py b/src/python_som/_enums.py index 2da5803..4c14925 100644 --- a/src/python_som/_enums.py +++ b/src/python_som/_enums.py @@ -1,25 +1,17 @@ """Names for the string-valued options, so a typo is a type error rather than a runtime one. -Every option these cover is still accepted as a plain string, and will be for the whole 0.4.x and -0.5.x series. ``mode=TrainingMode.BATCH`` and ``mode="batch"`` are interchangeable, compare equal, -hash equal, and serialise to the same JSON, because each member *is* a ``str``. - -**Both spellings are permanent.** 0.5.0 briefly deprecated plain strings and 0.6.0 withdrew that, -because every comparable library passes options as strings: scikit-learn -(``KMeans(init="k-means++")``), numpy (``np.pad(mode="constant")``), scipy -(``linkage(method="single")``), and both SOM peers, minisom and sompy. None of them export enums at -all. Being the only library in the ecosystem to reject ``mode="batch"`` would cost users more than -the consistency was worth. - -The enums remain because they cost nothing and some callers prefer them. The type-checking benefit -that motivated them is delivered by the ``Literal`` unions below rather than by removing anything: +Every option is also accepted as a plain string, permanently. ``mode=TrainingMode.BATCH`` and +``mode="batch"`` are interchangeable, compare equal, hash equal and serialise identically, because +each member *is* a ``str``. 0.5.0 briefly deprecated the string form and 0.6.0 withdrew that; the +changelog has the reasoning. + +The type-checking benefit comes from the ``Literal`` unions below rather than from the enums: ``mode="bacth"`` is a type error while ``mode="batch"`` is not. -**On the base class.** ``enum.StrEnum`` arrived in Python 3.11 and this package supports 3.10, so -:class:`_StrEnum` reproduces it. A bare ``class X(str, Enum)`` is *not* equivalent: its ``str()`` -returns ``'X.MEMBER'`` rather than the value, which would put the wrong text into any f-string, -filename or log line built from a member. Defining ``__str__`` explicitly makes the behaviour -identical on every supported version, which was checked on 3.10, 3.12 and 3.13 rather than assumed. +**On the base class.** ``enum.StrEnum`` needs Python 3.11 and this package supports 3.10, so +:class:`_StrEnum` reproduces it. A bare ``class X(str, Enum)`` is not equivalent: its ``str()`` +returns ``'X.MEMBER'`` rather than the value, which would put the wrong text into any f-string or +filename built from a member. """ from __future__ import annotations diff --git a/src/python_som/_som.py b/src/python_som/_som.py index 26886e6..cee3eb7 100644 --- a/src/python_som/_som.py +++ b/src/python_som/_som.py @@ -20,6 +20,7 @@ import numpy as np import numpy.typing as npt +from ._accelerate import bmu_kernel from ._artifact import ( ArtifactError, SOMConfig, @@ -36,12 +37,12 @@ from ._core._initialize import linear_models, random_models, sample_models from ._core._linalg import auto_dimensions from ._core._maps import activation_matrix, label_map, u_matrix, winner_map -from ._core._match import accumulate, activate, quantization, winner +from ._core._match import accumulate, activate, bmu_indices, quantization, winner from ._core._neighborhood import ( SIGNED_NEIGHBORHOODS, - kernel_view, + axis_matrix, resolve, - resolve_kernel, + resolve_axis_profile, ) from ._core._update import batch_update, stepwise_update from ._enums import ( @@ -146,19 +147,13 @@ def _warn_on_major_version_change(saved: str | None, path: object) -> None: def _validate_learning_rate(learning_rate: float) -> None: """Reject a learning rate that cannot train, and warn about one that is merely unwise. - Unchecked through 0.3.0, and the two failure modes are different in kind: - - A **non-positive** rate is rejected. ``alpha = 0`` freezes every model, so training runs to - completion and changes nothing. ``alpha = -1`` is worse: it moves models *away* from the samples - they match, taking the quantization error from 0.0 to 11.7 and the largest weight to 30 on a map - that started inside the unit cube. Neither can be what a caller meant, and both are silent. + A **non-positive** rate is rejected: ``alpha = 0`` freezes every model, and ``alpha = -1`` moves + them away from the samples they match, taking the quantization error from 0.0 to 11.7. Both are + silent failures. A rate **above 1** is warned about, not rejected. Eq. (3) moves a model a fraction ``alpha * h`` - of the way to the sample, so above 1 it overshoots and oscillates around the target rather than - settling on it. It does not necessarily diverge: measured at ``alpha = 5`` with decay disabled, - the largest weight stayed at 3.61, because the neighborhood damps the correction away from the - winner. Kohonen sets no upper bound, so rejecting it would invent a limit the sources do not - give. + of the way to the sample, so above 1 it overshoots and oscillates. It need not diverge, since + the neighborhood damps the correction away from the winner, and Kohonen gives no upper bound. :param learning_rate: The rate to check. :raises ValueError: If the rate is not a finite positive number. @@ -179,16 +174,6 @@ def _validate_learning_rate(learning_rate: float) -> None: class SOM: """A 2-D self-organizing map over NumPy arrays, pandas DataFrames or plain lists. - Features: - - Stepwise and batch training - - Random, random-sampling and linear (PCA) weight initialization - - Automatic selection of the map size ratio (with PCA) - - Support for cyclic arrays, for toroidal maps - - Gaussian, bubble and mexican hat neighborhood functions - - Support for custom decay functions - - Support for visualization (U-matrix, activation matrix) - - Support for supervised learning (label map) - Reference: Teuvo Kohonen, Essentials of the self-organizing map, Neural Networks 37 (2013) 52-65, https://doi.org/10.1016/j.neunet.2012.09.018 @@ -449,6 +434,10 @@ def train( "denominator is not sign-definite. Use mode='random' or mode='sequential'." ) raise ValueError(msg) + # A neighborhood that is unsigned but not separable would pass the check above and then be + # refused by `resolve_axis_profile` in `_train_batch`, which names the constraint. No + # separate guard here: every unsigned neighborhood is separable today, so it would be an + # untested branch for a case that cannot yet arise. if n_iteration is None: n_iteration = DEFAULT_ITERATIONS_PER_SAMPLE[mode] * len(array) @@ -518,13 +507,11 @@ def fit( ) -> SOM: """Train the map and return it, so calls can be chained. - ``y`` is accepted and ignored. Unsupervised estimators take it anyway, because that is what - lets ``Pipeline`` and ``cross_val_score`` call every step in the same way. + ``y`` is accepted and ignored, which is what lets ``Pipeline`` call every step the same way. - The training options are keyword arguments here rather than constructor arguments, so that - :class:`SOM` keeps one place where training is configured. The scikit-learn adapter in - :mod:`python_som.sklearn` takes them at construction instead, because ``get_params`` has to - expose them for ``GridSearchCV`` to tune them. + Training options are keyword arguments here rather than constructor arguments; + :mod:`python_som.sklearn` takes them at construction, because ``get_params`` must expose + them to ``GridSearchCV``. :param X: Training dataset of shape ``(n_samples, n_features)``. :param y: Ignored. @@ -567,20 +554,15 @@ def fit_transform( def predict(self, X: DataLike) -> npt.NDArray[np.integer]: # noqa: N803 """Return the index of the best-matching node for each sample. - A **flat** index, not a ``(row, column)`` pair. A 1-D array of labels is what scorers, - ``confusion_matrix`` and ``cross_val_score`` all assume, so returning coordinates would read - better for a grid and compose with nothing. Recover the grid position with - ``np.unravel_index(som.predict(X), som.get_shape())``, or call :meth:`winner` for a single - sample, which still returns ``(row, column)``. + A **flat** index, not a ``(row, column)`` pair, because that is what scorers and + ``confusion_matrix`` assume. Recover the grid position with + ``np.unravel_index(som.predict(X), som.get_shape())``; :meth:`winner` still returns + coordinates for a single sample. :param X: Dataset of shape ``(n_samples, n_features)``. :return: One flat node index per sample. """ - array = to_numpy(X) - shape = self._shape - return np.array( - [np.ravel_multi_index(self.winner(sample), shape) for sample in array], dtype=int - ) + return bmu_indices(to_numpy(X), self._weights, self._distance_function, bmu_kernel()) def score(self, X: DataLike, y: object = None) -> float: # noqa: ARG002, N803 """Return the negated quantization error, so that larger is better. @@ -623,12 +605,9 @@ def get_params(self, *, deep: bool = True) -> dict[str, Any]: # noqa: ARG002 def set_params(self, **params: Any) -> SOM: # noqa: ANN401 """Set constructor-level parameters in place and return this map. - Only the parameters that can be changed without rebuilding the models are accepted: the - rates, the radii and the decays. Changing the grid shape or ``input_len`` would invalidate - the weights, so those raise rather than silently leaving a map whose models do not match its - own description. - - This is what :meth:`set_learning_rate` and :meth:`set_neighborhood_radius` will become; both + Only what can change without rebuilding the models: the rates, radii and decays. The grid + shape and ``input_len`` raise, rather than leaving a map whose models do not match its own + description. Replaces :meth:`set_learning_rate` and :meth:`set_neighborhood_radius`, which still work and are removed in 1.0.0. :param params: Parameters to set. @@ -711,13 +690,11 @@ def config(self) -> SOMConfig: def save_npz(self, path: str | os.PathLike[str]) -> None: """Write the models and their provenance to a single ``.npz`` file. - The file holds the weights as an array and everything else as JSON beside them: the - configuration, the seed, the generator's current state, and the last training report. No - pickle is involved on either side, so the result is safe to load without executing code. + Weights as an array, everything else as JSON beside them: configuration, seed, generator + state and the last training report. No pickle on either side. - Saving the generator *state* as well as the seed is what lets :meth:`load_npz` resume the - same random stream. Re-seeding would restart it, and a resumed run would then silently - diverge from an uninterrupted one. + The generator *state* is saved as well as the seed, which is what lets :meth:`load_npz` + resume the same stream rather than restarting it. :param path: Destination file. """ @@ -748,13 +725,11 @@ def load_npz( ) -> SOM: """Rebuild a map saved by :meth:`save_npz`, models, generator state and all. - Continuing to train a loaded map produces the same weights as never having stopped, which is - the only useful definition of "loaded" for a stochastic process and is what the saved - generator state is for. + Continuing to train a loaded map gives the same weights as never having stopped, which is + what the saved generator state is for. - The four keyword arguments exist for maps trained with a function this package cannot look - up by name. Passing one the file did not need is harmless: it takes precedence over the - registered function of the same role. + The four keyword arguments are for maps trained with a function this package cannot resolve + by name. Passing one the file did not need is harmless. :param path: File to read. :param neighborhood_function: Replacement for a neighborhood that cannot be resolved. @@ -843,14 +818,8 @@ def _train_stepwise( ) -> tuple[float | None, float]: """Train one sample at a time, updating the winner and its neighbourhood. - Implements Eq. (3) of Kohonen (2013). ``'sequential'`` cycles through the dataset in order, - wrapping around until ``n_iteration`` steps have run. - - ``'random'`` draws samples **with replacement**, i.i.d., which is the stochastic - approximation of Robbins and Monro (1951) that Kohonen cites in Section 4.1. Before 0.3.0 - the draw used ``replace=(n_iteration > len(data))``, so it was a random permutation when the - iteration count did not exceed the sample count and i.i.d. only beyond it. That made the - character of the sampling depend on the iteration count, which is why it is now uniform. + Eq. (3) of Kohonen (2013). ``'sequential'`` cycles the dataset in order; ``'random'`` draws + i.i.d. **with replacement**, the Robbins-Monro approximation Kohonen cites in Section 4.1. :param array: Training dataset. :param n_iteration: Number of iterations. @@ -881,44 +850,26 @@ def _train_batch( ) -> tuple[float | None, float]: """Train with the batch algorithm, updating every model concurrently. - Implements Eq. (8) of Kohonen (2013). The winner map is recomputed from the models as they - stood at the start of each iteration, which is what makes the update concurrent. + Eq. (8) of Kohonen (2013). The winner map is recomputed from the models as they stood at + the start of each iteration, which is what makes the update concurrent (Section 4.4). - The neighborhood is evaluated **once per iteration**, not once per node. Eq. (8) needs - ``h_ji`` for every pair of nodes, and a neighborhood depends only on the offset between the - two -- so a single kernel over every offset serves the whole grid, and each node's - neighborhood is a slice of it. Evaluating per node instead made the neighborhood 42% of - batch training on a 40x40 map, more than the contraction it feeds; measured end to end, the - kernel is worth **1.2x to 1.5x**, more with the gaussian than the cheaper bubble. See - :func:`~python_som._core._neighborhood.offset_span` for why the offset-only dependence holds - on a torus as well as a flat grid, and ``benchmarks/bench_batch.py`` for the measurement. + The neighborhood is contracted as two per-axis matrices rather than evaluated per node; see + :doc:`/explanation/how-batch-training-is-computed`. :param array: Training dataset. :param n_iteration: Number of iterations. :param verbose: Whether to show a progress bar. """ - build_kernel = resolve_kernel(self._neighborhood_function_name) + profile = resolve_axis_profile(self._neighborhood_function_name) sigma = self._neighborhood_radius for t in self._progress(range(n_iteration), n_iteration, verbose=verbose): sigma = self._sigma(t, n_iteration) - sums, counts = accumulate(array, self._weights, self._shape, self._distance_function) - kernel = build_kernel(self._shape, sigma, self._cyclic) - - def neighborhood_of( - node: tuple[int, int], evaluated: npt.NDArray[np.floating] = kernel - ) -> npt.NDArray[np.floating]: - """Take this iteration's neighborhood for ``node`` out of the kernel. - - ``evaluated`` is a default argument rather than a closure over ``kernel`` so that - the value is bound at definition time, once per iteration. - - :param node: Coordinates of the node whose neighborhood is wanted. - :param evaluated: This iteration's kernel. - :return: Neighborhood weights over the grid, as a view into the kernel. - """ - return kernel_view(evaluated, self._shape, node) - - self._weights = batch_update(self._weights, sums, counts, neighborhood_of, self._shape) + sums, counts = accumulate( + array, self._weights, self._shape, self._distance_function, bmu_kernel() + ) + hx = axis_matrix(self._shape[0], sigma, cyclic=self._cyclic[0], profile=profile) + hy = axis_matrix(self._shape[1], sigma, cyclic=self._cyclic[1], profile=profile) + self._weights = batch_update(self._weights, sums, counts, hx, hy) # No learning rate: Eq. (8) is a weighted mean, so there is no step size to report. None # rather than the unused initial value, which would read as though it had been applied. diff --git a/src/python_som/_version.py b/src/python_som/_version.py index ba99836..d4ec072 100644 --- a/src/python_som/_version.py +++ b/src/python_som/_version.py @@ -9,4 +9,4 @@ __all__ = ["__version__"] -__version__ = "0.6.1" +__version__ = "0.7.0" diff --git a/src/python_som/sklearn.py b/src/python_som/sklearn.py index ffb0ec8..c00de49 100644 --- a/src/python_som/sklearn.py +++ b/src/python_som/sklearn.py @@ -1,24 +1,16 @@ """scikit-learn adapter. Import this only if you want a map to work inside scikit-learn. -:class:`~python_som.SOM` already provides ``fit``, ``transform``, ``predict`` and ``score``, which -is enough when *you* are the one calling them. It is not enough when scikit-learn does the calling: -since 1.7, ``Pipeline.predict``, ``GridSearchCV`` and ``cross_val_score`` all reach for -``__sklearn_tags__``, and the recommended way to have it is to inherit ``BaseEstimator``. +:class:`~python_som.SOM` already has ``fit``, ``transform``, ``predict`` and ``score``, which is +enough when you call them yourself. It is not enough when scikit-learn does: since 1.7, +``Pipeline.predict``, ``GridSearchCV`` and ``cross_val_score`` reach for ``__sklearn_tags__``, and +inheriting ``BaseEstimator`` is the supported way to have it. Defining that attribute by hand +couples to an internal that changed shape once already, and scikit-learn discourages it. -Measured against scikit-learn 1.9, the methods on :class:`~python_som.SOM` alone give ``clone`` and -``Pipeline.fit`` and then fail: ``Pipeline.predict``, ``GridSearchCV`` and ``cross_val_score`` all -raise ``AttributeError``. :class:`SOMEstimator` passes all five. - -So the integration lives here rather than in the core, and scikit-learn stays optional:: +So the integration lives here and scikit-learn stays optional:: pip install "python-som[sklearn]" -This is the ports-and-adapters shape the package already uses. An adapter may depend on the thing it -adapts; the core stays numpy-only, and importing :mod:`python_som` pulls none of this in. - -Defining ``__sklearn_tags__`` by hand was the alternative and was rejected: it couples to an -internal that already changed shape once between 1.6 and 1.7, and scikit-learn's own error message -says it does not recommend the approach. +The core stays numpy-only, and importing :mod:`python_som` pulls none of this in. """ from __future__ import annotations @@ -54,22 +46,18 @@ class SOMEstimator(ClusterMixin, TransformerMixin, BaseEstimator): """A self-organizing map as a scikit-learn estimator. - A SOM is a topologically-constrained k-means, so this follows ``KMeans``: ``transform`` gives a - cluster-distance space, ``predict`` gives one label per sample, ``score`` is negated so that - larger is better, and fitted attributes carry a trailing underscore. + Follows ``KMeans``: ``transform`` gives a cluster-distance space, ``predict`` one label per + sample, ``score`` a negated error so larger is better, fitted attributes a trailing underscore. >>> from python_som.sklearn import SOMEstimator >>> from sklearn.model_selection import GridSearchCV >>> search = GridSearchCV(SOMEstimator(), {"x": [4, 6]}, cv=3) # doctest: +SKIP - **Every argument is stored unmodified.** scikit-learn's ``clone`` rebuilds an estimator by - passing ``get_params()`` back to ``__init__`` and then checks the result is identical, so an - ``__init__`` that validates, coerces or derives anything breaks cloning. All of that is deferred - to :meth:`fit`, which is why this class holds settings rather than a - :class:`~python_som.SOM`. + **Every argument is stored unmodified.** ``clone`` rebuilds an estimator from ``get_params()`` + and checks the result is identical, so validating or deriving anything in ``__init__`` breaks + it. That is deferred to :meth:`fit`, which is why this class holds settings rather than a map. - ``input_len`` is deliberately absent: scikit-learn infers the feature count from ``X``, and - :attr:`n_features_in_` reports it after fitting. + ``input_len`` is absent: scikit-learn infers the feature count from ``X``. """ def __init__( @@ -94,8 +82,8 @@ def __init__( """Record the settings a map will be built from. Keyword-only after the grid dimensions, as ``KMeans`` is. The decays and the distance - default to None rather than to the functions themselves, so that the recorded parameters - stay exactly what the caller passed; :meth:`fit` substitutes the real defaults. + default to None so the recorded parameters stay exactly what the caller passed; + :meth:`fit` substitutes the real ones. :param x: Number of rows. :param y: Number of columns. @@ -133,10 +121,9 @@ def __init__( def fit(self, X: Any, y: Any = None, **kwargs: Any) -> SOMEstimator: # noqa: ANN401, ARG002, N803 """Build a map from the recorded settings and train it on ``X``. - A fresh map every call, unlike :meth:`python_som.SOM.fit`, which continues from wherever - its models were. Refitting an estimator is expected to start over: ``GridSearchCV`` fits the - same cloned estimator on fold after fold, and carrying weights between folds would leak one - fold into the next. + A fresh map every call, unlike :meth:`python_som.SOM.fit`, which continues. + ``GridSearchCV`` fits one cloned estimator fold after fold, and carrying weights over would + leak one fold into the next. :param X: Training dataset of shape ``(n_samples, n_features)``. :param y: Ignored. diff --git a/tests/test_batch_update_equivalence.py b/tests/test_batch_update_equivalence.py new file mode 100644 index 0000000..19d7102 --- /dev/null +++ b/tests/test_batch_update_equivalence.py @@ -0,0 +1,322 @@ +"""Batch training contracts the neighborhood by axis. It must equal evaluating it per node. + +Eq. (8) sums ``h`` over every pair of nodes. Because ``h`` depends only on the offset between two +nodes, that sum is a convolution, and a separable ``h`` turns it into two matrix contractions with +no loop over nodes. The saving is large, so the equality has to be held down hard rather than +assumed. + +Separability is an identity for exactly the two neighborhoods batch training admits, and for +neither of the reasons is it general: + +- the gaussian, because ``exp(-(dx^2 + dy^2) / 2s^2) == exp(-dx^2 / 2s^2) * exp(-dy^2 / 2s^2)``; +- the bubble, because ``max(|dx|, |dy|) <= r`` is the conjunction of two per-axis tests. + +The mexican hat factors under neither, which is why it has no axis profile and why batch training +rejects it. An outer product of two 1-D Ricker wavelets is a different function, positive in the +diagonal quadrants where the mexican hat must inhibit; that was a real defect in this package once, +and these tests are what stop the contraction quietly reintroducing it. +""" + +from __future__ import annotations + +import itertools + +import numpy as np +import pytest + +import python_som +from python_som import Neighborhood +from python_som._core._match import accumulate +from python_som._core._neighborhood import ( + AXIS_PROFILES, + NEIGHBORHOOD_FUNCTIONS, + SIGNED_NEIGHBORHOODS, + axis_matrix, + bubble, + gaussian, + resolve_axis_profile, +) +from python_som._core._update import batch_update + +#: Round-off scale for the contraction against the per-node reference. Measured at 3.1e-15 relative +#: on a 60x60 map; the two sum the same terms in a different order, so exact equality is not +#: available and asserting it would be asserting the wrong thing. +TOLERANCE = 1e-12 + +#: Fixed so a failure is reproducible. +SEED = 20260730 + +#: Shapes, including the degenerate single-row and single-column maps where an axis has length 1. +SHAPES = [(1, 6), (6, 1), (5, 5), (7, 4), (12, 9), (20, 16)] + +#: Radii, including one below 1 and one larger than the grid. +RADII = [0.5, 1.0, 2.5, 4.0, 30.0] + +#: Neighborhoods batch training admits. Derived rather than listed, so a new one joins the sweep. +SEPARABLE = sorted(AXIS_PROFILES) + + +def _per_node_reference( + weights: np.ndarray, + sums: np.ndarray, + counts: np.ndarray, + shape: tuple[int, int], + name: str, + sigma: float, + cyclic: tuple[bool, bool], +) -> np.ndarray: + """Evaluate Eq. (8) node by node from the isotropic definition. + + This is the definition the contraction has to match: it calls the public neighborhood function, + which is a function of ``sqdist``, once per node. + + :param weights: Current models. + :param sums: Per-node sums. + :param counts: Per-node counts. + :param shape: Grid shape. + :param name: Neighborhood function name. + :param sigma: Neighborhood radius. + :param cyclic: Whether each axis wraps. + :return: The updated models. + """ + evaluate = NEIGHBORHOOD_FUNCTIONS[name] + updated = weights.copy() + for node in np.ndindex(shape): + node_2d = (int(node[0]), int(node[1])) + h = evaluate(shape, node_2d, sigma, cyclic) + denominator = float(np.sum(h * counts)) + if denominator > 0: + updated[node_2d] = np.einsum("xy,xyf->f", h, sums) / denominator + return updated + + +def _case(shape: tuple[int, int], n_features: int = 3) -> tuple[np.ndarray, ...]: + """Build models, per-node sums and per-node counts for one grid. + + Counts are drawn with zeros in them on purpose: a node with no data in reach is the case that + must keep its previous value, and it is the one a naive implementation destroys. + + :param shape: Grid shape. + :param n_features: Number of features. + :return: Weights, sums and counts. + """ + rng = np.random.default_rng(SEED + shape[0] * 100 + shape[1]) + return ( + rng.normal(size=(*shape, n_features)), + rng.normal(size=(*shape, n_features)), + rng.integers(0, 3, size=shape).astype(float), + ) + + +# --------------------------------------------------------------------------------------------- +# The contraction equals the definition +# --------------------------------------------------------------------------------------------- + + +@pytest.mark.parametrize("shape", SHAPES) +@pytest.mark.parametrize("name", SEPARABLE) +@pytest.mark.parametrize("cyclic", list(itertools.product([False, True], repeat=2))) +@pytest.mark.parametrize("sigma", RADII) +def test_the_contraction_equals_the_per_node_definition( + shape: tuple[int, int], name: str, cyclic: tuple[bool, bool], sigma: float +) -> None: + """Every shape, every neighborhood, every cyclic combination, every radius.""" + weights, sums, counts = _case(shape) + profile = resolve_axis_profile(name) + hx = axis_matrix(shape[0], sigma, cyclic=cyclic[0], profile=profile) + hy = axis_matrix(shape[1], sigma, cyclic=cyclic[1], profile=profile) + + contracted = batch_update(weights, sums, counts, hx, hy) + reference = _per_node_reference(weights, sums, counts, shape, name, sigma, cyclic) + + scale = max(float(np.abs(reference).max()), 1.0) + assert float(np.abs(contracted - reference).max()) / scale < TOLERANCE + + +@pytest.mark.parametrize("name", SEPARABLE) +@pytest.mark.parametrize("cyclic", list(itertools.product([False, True], repeat=2))) +def test_the_axis_factors_multiply_to_the_isotropic_neighborhood( + name: str, cyclic: tuple[bool, bool] +) -> None: + """The claim separability rests on, asserted directly rather than only through Eq. (8). + + For each node, the outer product of the two axis factors must be that node's neighborhood as the + public function computes it from ``sqdist``. + """ + shape, sigma = (9, 7), 2.0 + profile = resolve_axis_profile(name) + hx = axis_matrix(shape[0], sigma, cyclic=cyclic[0], profile=profile) + hy = axis_matrix(shape[1], sigma, cyclic=cyclic[1], profile=profile) + evaluate = NEIGHBORHOOD_FUNCTIONS[name] + + for node in np.ndindex(shape): + node_2d = (int(node[0]), int(node[1])) + factored = np.multiply.outer(hx[:, node_2d[0]], hy[:, node_2d[1]]) + np.testing.assert_allclose(factored, evaluate(shape, node_2d, sigma, cyclic), atol=1e-15) + + +def test_a_node_with_no_data_in_reach_keeps_its_previous_value() -> None: + """Kohonen Eq. (8) is undefined where the denominator is zero, so the old model stands. + + Regression for a defect that wiped 282 of 900 models in a single step on a 30x30 map by building + the result from a zeroed array. The bubble makes it reachable: it is exactly zero outside its + radius, where the gaussian is merely small. + """ + shape, sigma = (30, 30), 1.0 + weights, sums, counts = _case(shape) + counts[:] = 0.0 + counts[0, 0] = 5.0 + + profile = resolve_axis_profile("bubble") + hx = axis_matrix(shape[0], sigma, cyclic=False, profile=profile) + hy = axis_matrix(shape[1], sigma, cyclic=False, profile=profile) + updated = batch_update(weights, sums, counts, hx, hy) + + reached = np.zeros(shape, dtype=bool) + reached[:2, :2] = True + np.testing.assert_array_equal(updated[~reached], weights[~reached]) + assert not np.array_equal(updated[0, 0], weights[0, 0]), "the node with data must have moved" + + +def test_the_update_is_concurrent_over_every_node() -> None: + """Kohonen Section 4.4: models are replaced "in one concurrent computing operation". + + Every node must be computed from the models as they stood at the start of the iteration. A loop + writing into the array it reads would satisfy the other tests here and fail this one. + """ + shape, sigma = (8, 6), 2.0 + weights, sums, counts = _case(shape) + profile = resolve_axis_profile("gaussian") + hx = axis_matrix(shape[0], sigma, cyclic=False, profile=profile) + hy = axis_matrix(shape[1], sigma, cyclic=False, profile=profile) + + updated = batch_update(weights, sums, counts, hx, hy) + + # Recompute one late node from the *original* models. If anything had leaked from an earlier + # node's new value, this would disagree. + late = (shape[0] - 1, shape[1] - 1) + h = gaussian(shape, late, sigma, (False, False)) + expected = np.einsum("xy,xyf->f", h, sums) / float(np.sum(h * counts)) + np.testing.assert_allclose(updated[late], expected, rtol=1e-12) + + assert not np.shares_memory(updated, weights), "the update must not alias its input" + + +# --------------------------------------------------------------------------------------------- +# The registry is the batch-legality rule +# --------------------------------------------------------------------------------------------- + + +def test_only_separable_neighborhoods_have_an_axis_profile() -> None: + """The mexican hat must never acquire one. + + ``(1 - u) exp(-u)`` does not factor. An outer product of two 1-D Ricker wavelets is a different + function: it is positive in the diagonal quadrants, +0.165 at 2 sigma where the correct value is + -0.055, placing an excitatory lobe where the mexican hat must inhibit. + """ + assert set(AXIS_PROFILES) == {"gaussian", "bubble"} + assert SIGNED_NEIGHBORHOODS.isdisjoint(AXIS_PROFILES) + + +def test_every_unsigned_neighborhood_is_separable() -> None: + """What batch training relies on: anything it accepts, the contraction can express. + + If a future neighborhood is unsigned but not separable, this fails and the choice becomes + explicit rather than silently approximated. + """ + unsigned = set(NEIGHBORHOOD_FUNCTIONS) - set(SIGNED_NEIGHBORHOODS) + assert unsigned == set(AXIS_PROFILES) + + +def test_resolve_axis_profile_rejects_a_non_separable_neighborhood() -> None: + """The mexican hat reaches this only if the signed check goes; the message still names why.""" + with pytest.raises(ValueError, match="not separable"): + resolve_axis_profile("mexican_hat") + + +def test_resolve_axis_profile_rejects_an_unknown_name() -> None: + with pytest.raises(ValueError, match="not separable"): + resolve_axis_profile("spectral") + + +@pytest.mark.parametrize("name", SEPARABLE) +def test_the_axis_profiles_validate_the_radius(name: str) -> None: + """Validation lives in the profile, so the contraction path cannot skip it.""" + profile = resolve_axis_profile(name) + with pytest.raises(ValueError, match="must be a finite"): + profile(np.array([0.0, 1.0]), float("nan")) + with pytest.raises(ValueError, match="must be a finite"): + profile(np.array([0.0, 1.0]), -1.0) + + +def test_the_bubble_accepts_a_zero_radius_and_the_gaussian_does_not() -> None: + """Unchanged from the per-node forms: zero selects the winner alone, or divides by zero.""" + np.testing.assert_array_equal( + AXIS_PROFILES["bubble"](np.array([-1.0, 0.0, 1.0]), 0.0), np.array([0.0, 1.0, 0.0]) + ) + with pytest.raises(ValueError, match="must be a finite positive"): + AXIS_PROFILES["gaussian"](np.array([0.0]), 0.0) + + +# --------------------------------------------------------------------------------------------- +# End to end +# --------------------------------------------------------------------------------------------- + + +@pytest.mark.parametrize("neighborhood", [Neighborhood.GAUSSIAN, Neighborhood.BUBBLE]) +@pytest.mark.parametrize("cyclic", list(itertools.product([False, True], repeat=2))) +def test_batch_training_matches_the_per_node_definition_end_to_end( + neighborhood: Neighborhood, cyclic: tuple[bool, bool] +) -> None: + """A whole training run, not one update, so any per-iteration drift accumulates into view.""" + shape, n_iteration = (12, 9), 20 + rng = np.random.default_rng(SEED) + data = rng.normal(size=(90, 4)) + initial = rng.normal(size=(*shape, 4)) + + som = python_som.SOM( + x=shape[0], + y=shape[1], + input_len=4, + neighborhood_function=neighborhood, + neighborhood_radius=3.0, + cyclic_x=cyclic[0], + cyclic_y=cyclic[1], + random_seed=SEED, + ) + som._weights = initial.copy() + som.train(data, n_iteration=n_iteration, mode="batch") + + reference = initial.copy() + for step in range(n_iteration): + sigma = som._sigma(step, n_iteration) + sums, counts = accumulate(data, reference, shape, som._distance_function) + reference = _per_node_reference( + reference, sums, counts, shape, neighborhood.value, sigma, cyclic + ) + + scale = float(np.abs(reference).max()) + assert float(np.abs(som.get_weights() - reference).max()) / scale < TOLERANCE + + +def test_batch_training_still_rejects_the_mexican_hat() -> None: + """Unchanged, and the message is still about the sign rather than about separability.""" + som = python_som.SOM(x=6, y=6, input_len=3, neighborhood_function="mexican_hat", random_seed=1) + with pytest.raises(ValueError, match="cannot be used with the 'batch' training mode"): + som.train(np.random.default_rng(0).normal(size=(20, 3)), n_iteration=5, mode="batch") + + +def test_the_bubble_is_not_isotropic_under_the_euclidean_metric() -> None: + """Its metric is Chebyshev, which is why it factors. Pinned because the sources disagree. + + Kohonen Section 4.2 describes the flat neighborhood as "1 up to a certain radius from the + winner", which reads Euclidean; Vrieze's appendix computes ``MAX(ABS(i - w_i), ABS(j - w_j))``, + which is Chebyshev, and that is what this package implements. A Euclidean disc would not be + separable and so could not use this path at all. + + The smallest counterexample is a radius of ``sqrt(50)``: ``(5, 5)`` lies inside a ``sigma = 5`` + square while ``(7, 1)`` lies outside, at equal Euclidean distance from the winner. + """ + h = bubble((15, 15), (0, 0), 5.0, (False, False)) + assert h[5, 5] == 1.0 + assert h[7, 1] == 0.0 diff --git a/tests/test_bmu_search.py b/tests/test_bmu_search.py new file mode 100644 index 0000000..e97b576 --- /dev/null +++ b/tests/test_bmu_search.py @@ -0,0 +1,253 @@ +"""The vectorised best-matching-unit search must select the same nodes as the definition. + +Eq. (4) of Kohonen (2013) is ``c = argmin_i ||x - m_i||``. For the Euclidean distance the search +expands that norm into a matrix product, which is much faster and is *not* obviously the same thing. +These tests hold the two properties that make it the same thing: it picks the node the definition +picks, and it stays the Euclidean map rather than becoming the dot-product map of Section 4.5. + +The expansion also has a failure mode that only appears far from the origin, and it is severe enough +to have its own regression test below. +""" + +from __future__ import annotations + +import numpy as np +import pytest + +import python_som +from python_som._core._distance import euclidean_distance +from python_som._core._match import accumulate, bmu_indices, quantization, winner + +#: Fixed so a failure is reproducible. +SEED = 20260730 + + +def _exact(data: np.ndarray, weights: np.ndarray) -> np.ndarray: + """Select the best-matching node by the definition, one full norm per sample. + + :param data: Dataset. + :param weights: Models. + :return: One flat node index per sample. + """ + flat = weights.reshape(-1, weights.shape[-1]) + return np.array([np.linalg.norm(x - flat, axis=-1).argmin() for x in data]) + + +def _case( + shape: tuple[int, int], n_samples: int, n_features: int, offset: float = 0.0 +) -> tuple[np.ndarray, np.ndarray]: + """Build models and a dataset, optionally far from the origin. + + :param shape: Grid shape. + :param n_samples: Number of samples. + :param n_features: Number of features. + :param offset: Constant added to both, to move them away from the origin. + :return: Models and dataset. + """ + rng = np.random.default_rng(SEED) + return ( + rng.normal(size=(*shape, n_features)) + offset, + rng.normal(size=(n_samples, n_features)) + offset, + ) + + +# --------------------------------------------------------------------------------------------- +# It selects what the definition selects +# --------------------------------------------------------------------------------------------- + + +@pytest.mark.parametrize( + ("shape", "n_samples", "n_features"), + [((5, 5), 50, 3), ((20, 20), 200, 4), ((40, 30), 500, 8), ((1, 9), 40, 2), ((60, 60), 300, 12)], +) +def test_the_fast_search_selects_the_same_nodes( + shape: tuple[int, int], n_samples: int, n_features: int +) -> None: + """Identical indices, not close ones: a different node is a different answer.""" + weights, data = _case(shape, n_samples, n_features) + np.testing.assert_array_equal( + bmu_indices(data, weights, euclidean_distance), _exact(data, weights) + ) + + +def test_it_agrees_with_winner_for_every_sample() -> None: + """The single-sample path and the whole-dataset path must not drift apart.""" + weights, data = _case((7, 5), 60, 4) + flat = bmu_indices(data, weights, euclidean_distance) + rows, columns = np.unravel_index(flat, weights.shape[:2]) + for sample, row, column in zip(data, rows, columns, strict=True): + assert (int(row), int(column)) == winner(sample, weights, euclidean_distance) + + +def test_ties_go_to_the_first_node_in_c_order() -> None: + """Arbitrary but fixed, and it must match ``argmin``, which is what ``winner`` uses.""" + weights = np.full((3, 3, 2), 10.0) + weights[0, 2] = [1.0, 0.0] + weights[2, 0] = [1.0, 0.0] + data = np.array([[1.0, 0.0]]) + + assert int(bmu_indices(data, weights, euclidean_distance)[0]) == 2 + assert winner(data[0], weights, euclidean_distance) == (0, 2) + + +def test_a_chunk_boundary_does_not_change_the_result() -> None: + """The search runs in blocks, so a dataset larger than one block exercises the seam. + + At 60x60 the block holds about 17 samples, so 500 samples cross it many times. + """ + weights, data = _case((60, 60), 500, 6) + np.testing.assert_array_equal( + bmu_indices(data, weights, euclidean_distance), _exact(data, weights) + ) + + +# --------------------------------------------------------------------------------------------- +# The cancellation the expansion would otherwise cause +# --------------------------------------------------------------------------------------------- + + +@pytest.mark.parametrize("offset", [0.0, 1e3, 1e6, 1e9, 1e12]) +def test_the_search_is_exact_far_from_the_origin(offset: float) -> None: + """Regression for catastrophic cancellation, with the measured numbers. + + ``||x - w||^2 = ||x||^2 - 2 x.w + ||w||^2`` is exact in real arithmetic and not in floating + point. With models offset by 1e9, ``||w||^2`` is about 1e18 while the differences between models + are of order 1, so the subtraction loses every significant digit. Measured without centring: + + ====== ========================= + offset samples given a wrong node + ====== ========================= + 1e6 0 of 500 + 1e9 **500 of 500** + 1e12 **500 of 500** + ====== ========================= + + Subtracting the models' mean from both sides is exact in ``||x - w||``, costs 1%, and removes it + at every offset above. Data far from the origin is not exotic: timestamps, easting and northing + coordinates and absolute sensor readings all look like this, and it is the same failure this + package fixed in linear initialization in 0.4.0. + """ + weights, data = _case((40, 40), 500, 6, offset=offset) + np.testing.assert_array_equal( + bmu_indices(data, weights, euclidean_distance), _exact(data, weights) + ) + + +def test_the_uncentred_expansion_really_does_fail_there() -> None: + """Show the defect the centring prevents, so the fix is not mistaken for a redundant line. + + Without this, someone reading ``flat - shift`` sees an operation with no visible effect and + deletes it, and every test above still passes at the offsets they happen to try. + """ + weights, data = _case((40, 40), 500, 6, offset=1e9) + flat = weights.reshape(-1, weights.shape[-1]) + + uncentred = (np.einsum("nf,nf->n", flat, flat)[None, :] - 2.0 * (data @ flat.T)).argmin(axis=1) + + wrong = int((uncentred != _exact(data, weights)).sum()) + assert wrong > len(data) // 2, ( + f"expected the uncentred expansion to fail badly at 1e9, got {wrong} of {len(data)}" + ) + + +# --------------------------------------------------------------------------------------------- +# It is still the Euclidean map +# --------------------------------------------------------------------------------------------- + + +def test_it_is_not_the_dot_product_map_of_section_4_5() -> None: + """Kohonen Section 4.5 defines a different algorithm, and this is not it. + + Eq. (9), ``c = argmax_i dot(x, m_i)``, requires the models to be "kept normalized to constant + length all the time" and selects a different node when they are not. A matrix product in the + winner search reads exactly like a silent switch to it, so the difference is asserted. + + The models here have deliberately unequal lengths, which is what makes the two criteria diverge. + """ + weights = np.array([[[1.0, 0.0], [10.0, 10.0]]]) + data = np.array([[1.0, 0.5]]) + + flat = weights.reshape(-1, 2) + assert int(np.argmax(flat @ data[0])) == 1, "the dot-product criterion prefers the long model" + assert int(bmu_indices(data, weights, euclidean_distance)[0]) == 0 + assert winner(data[0], weights, euclidean_distance) == (0, 0) + + +# --------------------------------------------------------------------------------------------- +# A custom distance keeps the exact path +# --------------------------------------------------------------------------------------------- + + +def _manhattan(x: object, weights: object) -> np.ndarray: + """Sum of absolute differences along the last axis. + + :param x: Input vector. + :param weights: One model or an array of them. + :return: Distances. + """ + result: np.ndarray = np.abs(np.asarray(x) - np.asarray(weights)).sum(axis=-1) + return result + + +def test_a_custom_distance_is_used_rather_than_the_fast_path() -> None: + """The expansion is an identity for the Euclidean norm alone, so anything else takes the loop. + + Constructed so the two metrics disagree: under Manhattan the first model wins, under Euclidean + the second does. If the fast path were taken regardless, this would return the Euclidean answer. + """ + weights = np.array([[[0.9, 0.9], [0.0, 1.4]]]) + data = np.array([[0.0, 0.0]]) + + assert int(bmu_indices(data, weights, _manhattan)[0]) == 1 + assert int(bmu_indices(data, weights, euclidean_distance)[0]) == 0 + + +def test_both_paths_agree_when_the_custom_distance_is_euclidean() -> None: + """A user-supplied function that happens to be Euclidean must give the same nodes. + + ``python_som.euclidean_distance`` is selected by identity, so passing an equivalent but distinct + callable takes the slow path. The two must still agree. + """ + weights, data = _case((10, 8), 120, 5) + + def same_but_not_identical(x: object, w: object) -> np.ndarray: + """Euclidean distance, written out so it is not the registered function object.""" + result: np.ndarray = np.linalg.norm(np.asarray(x) - np.asarray(w), axis=-1) + return result + + np.testing.assert_array_equal( + bmu_indices(data, weights, same_but_not_identical), + bmu_indices(data, weights, euclidean_distance), + ) + + +def test_accumulate_and_quantization_go_through_the_same_search() -> None: + """Eq. (8)'s inputs and the reported error must agree with the nodes the search chose.""" + shape = (9, 7) + weights, data = _case(shape, 150, 4) + + _, counts = accumulate(data, weights, shape, euclidean_distance) + nodes = bmu_indices(data, weights, euclidean_distance) + + expected_counts = np.bincount(nodes, minlength=shape[0] * shape[1]).reshape(shape) + np.testing.assert_array_equal(counts, expected_counts.astype(float)) + assert counts.sum() == len(data) + + flat = weights.reshape(-1, weights.shape[-1]) + errors = quantization(data, weights, euclidean_distance) + for error, sample, node in zip(errors, data, nodes, strict=True): + assert error == pytest.approx(float(np.linalg.norm(sample - flat[node]))) + + +def test_quantization_error_is_unchanged_by_the_faster_search() -> None: + """The reported number is a distance, not the search's score, which drops a constant term.""" + som = python_som.SOM(x=8, y=6, input_len=4, random_seed=SEED) + rng = np.random.default_rng(SEED) + data = rng.normal(size=(80, 4)) + som.weight_initialization(mode="random") + + flat = som.get_weights().reshape(-1, 4) + expected = float( + np.mean([np.linalg.norm(x - flat, axis=-1).min() for x in data]), + ) + assert som.quantization_error(data) == pytest.approx(expected) diff --git a/tests/test_core_boundary.py b/tests/test_core_boundary.py index d768a6b..4d7f1b5 100644 --- a/tests/test_core_boundary.py +++ b/tests/test_core_boundary.py @@ -17,7 +17,7 @@ import python_som from python_som import WeightInit from python_som._core import _update -from python_som._core._neighborhood import gaussian +from python_som._core._neighborhood import axis_matrix, gaussian, gaussian_axis_profile from tests.conftest import MODEL_SEED, make_som #: The core package on disk, scanned rather than imported. @@ -154,9 +154,9 @@ def test_batch_update_leaves_unreached_models_untouched() -> None: weights = np.arange(4 * 4 * 2, dtype=float).reshape((*shape, 2)) sums = np.zeros((*shape, 2)) counts = np.zeros(shape) # no data anywhere - result = _update.batch_update( - weights, sums, counts, lambda node: gaussian(shape, node, 1.0, (False, False)), shape - ) + hx = axis_matrix(shape[0], 1.0, cyclic=False, profile=gaussian_axis_profile) + hy = axis_matrix(shape[1], 1.0, cyclic=False, profile=gaussian_axis_profile) + result = _update.batch_update(weights, sums, counts, hx, hy) np.testing.assert_array_equal(result, weights) @@ -167,9 +167,9 @@ def test_batch_update_returns_a_new_array() -> None: original = weights.copy() counts = np.ones(shape) sums = np.full((*shape, 2), 5.0) - result = _update.batch_update( - weights, sums, counts, lambda node: gaussian(shape, node, 1.0, (False, False)), shape - ) + hx = axis_matrix(shape[0], 1.0, cyclic=False, profile=gaussian_axis_profile) + hy = axis_matrix(shape[1], 1.0, cyclic=False, profile=gaussian_axis_profile) + result = _update.batch_update(weights, sums, counts, hx, hy) assert result is not weights np.testing.assert_array_equal(weights, original) diff --git a/tests/test_kernel_equivalence.py b/tests/test_kernel_equivalence.py deleted file mode 100644 index b517350..0000000 --- a/tests/test_kernel_equivalence.py +++ /dev/null @@ -1,311 +0,0 @@ -"""The kernel form of each neighborhood must equal the per-node form exactly, not approximately. - -Batch training evaluates the neighborhood once per iteration and slices it per node, rather than -evaluating it once per node. That is only sound because a neighborhood depends on the offset between -two nodes and never on where the winner sits, so these tests exist to hold that property down. - -The bar is **exactly 0.0**, not a tolerance. A speedup that moves trained weights is a bug, and a -tolerance would hide precisely the class of error this replaces: two implementations of one formula -that agree on the cases someone thought to check. - -The design makes the equality structural rather than hoped for -- both forms call the same private -profile, differing only in which offsets they pass -- so these tests guard the *premise* -(offset-only dependence, and the right slice) rather than a typo in a second copy of a formula. -""" - -from __future__ import annotations - -import functools -import itertools -from typing import TYPE_CHECKING - -import numpy as np -import pytest - -import python_som -from python_som import Neighborhood -from python_som._core._match import accumulate -from python_som._core._neighborhood import ( - NEIGHBORHOOD_FUNCTIONS, - NEIGHBORHOOD_KERNELS, - axis_offsets, - bubble, - gaussian, - kernel_view, - mexican_hat, - offset_span, - resolve_kernel, -) -from python_som._core._update import batch_update - -if TYPE_CHECKING: # pragma: no cover - from collections.abc import Callable - - from python_som._core._neighborhood import NeighborhoodFunction - -#: Grid shapes to sweep, including degenerate single-row and single-column maps, where the offset -#: span collapses and an off-by-one in the slice would be invisible on a square grid. -SHAPES = [(10, 10), (7, 13), (20, 20), (9, 4), (1, 5), (6, 1)] - -#: Radii to sweep. ``0.0`` is admissible for the bubble alone, and is included because it is the one -#: value where the neighborhood is a single node and the slice has to be exactly right. -RADII = [0.0, 0.5, 1.0, 2.5, 4.0, 7.0] - -#: All four combinations, so a mixed toroidal map (one axis wrapping, one not) is covered. Each axis -#: folds independently, which is why one slice serves every combination. -CYCLIC = list(itertools.product([False, True], repeat=2)) - -#: The three distinct functions, ignoring the ``mexican_hat``/``mexicanhat`` alias. -FUNCTIONS = {"gaussian": gaussian, "bubble": bubble, "mexican_hat": mexican_hat} - - -def _evaluate_per_node( - function: NeighborhoodFunction, - shape: tuple[int, int], - sigma: float, - cyclic: tuple[bool, bool], - node: tuple[int, int], -) -> np.ndarray: - """Evaluate one node's neighborhood directly, as batch training did before the kernel. - - Takes everything explicitly at module level rather than closing over the loop variables, so a - late-binding mistake cannot quietly make both arms of the comparison the same. - - :param function: The per-node neighborhood function. - :param shape: Shape of the grid. - :param sigma: This iteration's radius. - :param cyclic: Whether each axis wraps. - :param node: Node whose neighborhood is wanted. - :return: Neighborhood weights over the grid. - """ - return function(shape, node, sigma, cyclic) - - -def _admissible(name: str, sigma: float) -> bool: - """Whether this function accepts this radius. - - :param name: Neighborhood function name. - :param sigma: Radius. - :return: True if the call would not raise. - """ - return sigma > 0 or name == "bubble" - - -@pytest.mark.parametrize("shape", SHAPES) -@pytest.mark.parametrize("cyclic", CYCLIC) -def test_kernel_equals_per_node_evaluation_for_every_node( - shape: tuple[int, int], cyclic: tuple[bool, bool] -) -> None: - """Sweep every function, radius and **node** of the grid, asserting exact equality. - - This is the load-bearing test of the optimization. Across all parameters it covers 40,832 - (function, shape, cyclic, radius, node) combinations, every one of which must agree at 0.0. - """ - for name, function in FUNCTIONS.items(): - build = resolve_kernel(name) - for sigma in RADII: - if not _admissible(name, sigma): - continue - kernel = build(shape, sigma, cyclic) - assert kernel.shape == (2 * shape[0] - 1, 2 * shape[1] - 1) - for node in itertools.product(range(shape[0]), range(shape[1])): - expected = function(shape, node, sigma, cyclic) - actual = kernel_view(kernel, shape, node) - difference = np.abs(expected - actual).max() - assert difference == 0.0, ( - f"{name} on {shape}, cyclic={cyclic}, sigma={sigma}, node={node}: " - f"kernel and per-node evaluation differ by {difference}" - ) - - -def test_the_sweep_really_covers_every_node() -> None: - """Guard the test above against silently shrinking. - - A parametrised sweep that stops covering what its docstring claims is worse than no sweep, so - the count is asserted rather than described. - """ - total = sum( - shape[0] * shape[1] - for shape in SHAPES - for cyclic in CYCLIC - for name in FUNCTIONS - for sigma in RADII - if _admissible(name, sigma) - ) - assert total == 40832, total - - -@pytest.mark.parametrize("name", sorted(NEIGHBORHOOD_KERNELS)) -def test_every_registered_function_has_a_kernel(name: str) -> None: - """Batch training takes the kernel path unconditionally, with no fallback branch. - - That is only safe if the two registries agree, so it is asserted rather than assumed. A name in - ``NEIGHBORHOOD_FUNCTIONS`` without a kernel would be an ``AttributeError`` deep in training. - """ - assert name in NEIGHBORHOOD_FUNCTIONS - assert callable(NEIGHBORHOOD_KERNELS[name]) - - -def test_the_two_registries_have_the_same_keys() -> None: - """Including the ``mexican_hat``/``mexicanhat`` alias, which is easy to add to only one.""" - assert set(NEIGHBORHOOD_KERNELS) == set(NEIGHBORHOOD_FUNCTIONS) - - -def test_kernel_view_is_a_view_and_not_a_copy() -> None: - """Copying the slice would give back much of what evaluating once saved. - - ``batch_update`` calls this for every node, so an accidental copy would allocate ``x * y`` - floats per node per iteration -- exactly the cost the kernel exists to avoid. - """ - kernel = NEIGHBORHOOD_KERNELS["gaussian"]((12, 9), 2.0, (False, False)) - view = kernel_view(kernel, (12, 9), (5, 4)) - assert np.shares_memory(view, kernel), "kernel_view must not copy" - assert view.base is not None - - -@pytest.mark.parametrize("cyclic", [False, True]) -def test_offset_span_covers_exactly_the_reachable_offsets(cyclic: bool) -> None: - """The span must contain every offset ``i - c`` that any pair of nodes can produce. - - One element short and the slice for a corner node would read out of bounds or silently wrap. - """ - length = 9 - span = offset_span(length, cyclic=cyclic) - assert span.shape == (2 * length - 1,) - - reachable = { - float(offset) - for centre in range(length) - for offset in axis_offsets(length, centre, cyclic=cyclic) - } - assert reachable <= set(span.tolist()) - - -@pytest.mark.parametrize("cyclic", [False, True]) -def test_offset_span_agrees_with_axis_offsets_elementwise(cyclic: bool) -> None: - """The span is ``axis_offsets`` read at a shifted origin, which is what makes the slice valid. - - Asserted per element so a fold applied with the wrong period would fail here rather than as a - puzzling difference in trained weights. This is the specific trap: the span is ``2L-1`` wide but - must fold with period ``L``. - """ - length = 11 - span = offset_span(length, cyclic=cyclic) - for centre in range(length): - expected = axis_offsets(length, centre, cyclic=cyclic) - actual = span[length - 1 - centre : 2 * length - 1 - centre] - np.testing.assert_array_equal(actual, expected) - - -def test_a_cyclic_span_cannot_be_faked_with_a_wider_axis() -> None: - """Pin the reason ``offset_span`` exists instead of reusing ``axis_offsets`` on a ``2L-1`` axis. - - On a flat grid the two coincide, which is what makes this an easy and wrong simplification: the - fold has to use the real period ``L``, and an axis of width ``2L-1`` folds with the wrong one. - """ - length = 10 - correct = offset_span(length, cyclic=True) - naive = axis_offsets(2 * length - 1, length - 1, cyclic=True) - assert not np.array_equal(correct, naive), ( - "if these ever agree, the simplification is safe and this test should be revisited" - ) - - -@pytest.mark.parametrize("sigma", [0.0, -1.0, -0.5, float("nan"), float("inf")]) -def test_kernels_validate_the_radius_exactly_as_the_per_node_form_does(sigma: float) -> None: - """Validation lives in the shared profile, so neither form can accept what the other rejects.""" - shape = (6, 6) - for name, function in FUNCTIONS.items(): - build = resolve_kernel(name) - per_node_raised = kernel_raised = False - try: - function(shape, (3, 3), sigma, (False, False)) - except ValueError: - per_node_raised = True - try: - build(shape, sigma, (False, False)) - except ValueError: - kernel_raised = True - assert per_node_raised == kernel_raised, ( - f"{name} at sigma={sigma}: per-node raised={per_node_raised}, kernel={kernel_raised}" - ) - - -def test_resolve_kernel_rejects_an_unknown_name() -> None: - """The same error shape as ``resolve``, naming the valid options.""" - with pytest.raises(ValueError, match="Invalid value for 'neighborhood_function' parameter"): - resolve_kernel("spectral") - - -# --------------------------------------------------------------------------------------------- -# End to end: training through the kernel equals training through per-node evaluation -# --------------------------------------------------------------------------------------------- - - -@pytest.mark.parametrize("neighborhood", [Neighborhood.GAUSSIAN, Neighborhood.BUBBLE]) -@pytest.mark.parametrize("cyclic", [(False, False), (True, True), (True, False)]) -def test_batch_training_is_unchanged_by_the_kernel( - neighborhood: Neighborhood, cyclic: tuple[bool, bool] -) -> None: - """Many iterations, with a decaying radius, against the per-node path it replaced. - - The single-iteration case is already covered by - ``test_batch_matches_a_reference_implementation`` in ``tests/test_training.py``, which builds - Eq. (8) as a literal double loop. This one runs the - full loop for long enough that the radius decays through several values, because the kernel is - rebuilt on each iteration and a stale-kernel bug would only show up after the first. - - Only gaussian and bubble appear: batch training rejects signed neighborhoods, so the mexican hat - cannot reach this path at all. - """ - rng = np.random.default_rng(20260731) - data = rng.normal(size=(80, 3)) - - weights = np.asarray( - python_som.SOM( - x=7, - y=5, - input_len=3, - neighborhood_function=neighborhood, - neighborhood_radius=3.0, - cyclic_x=cyclic[0], - cyclic_y=cyclic[1], - random_seed=11, - ).get_weights() - ) - - def train(*, use_kernel: bool) -> np.ndarray: - """Run batch training with either the kernel path or a per-node one. - - :param use_kernel: Whether to slice a kernel or evaluate the neighborhood per node. - :return: The trained models. - """ - som = python_som.SOM( - x=7, - y=5, - input_len=3, - neighborhood_function=neighborhood, - neighborhood_radius=3.0, - cyclic_x=cyclic[0], - cyclic_y=cyclic[1], - random_seed=11, - ) - som._weights = weights.copy() - shape = som.get_shape() - current = weights.copy() - for step in range(12): - sigma = som._sigma(step, 12) - sums, counts = accumulate(data, current, shape, som._distance_function) - neighborhood_of: Callable[[tuple[int, int]], np.ndarray] - if use_kernel: - kernel = resolve_kernel(neighborhood.value)(shape, sigma, cyclic) - neighborhood_of = functools.partial(kernel_view, kernel, shape) - else: - neighborhood_of = functools.partial( - _evaluate_per_node, FUNCTIONS[neighborhood.value], shape, sigma, cyclic - ) - current = batch_update(current, sums, counts, neighborhood_of, shape) - return current - - difference = np.abs(train(use_kernel=True) - train(use_kernel=False)).max() - assert difference == 0.0, f"kernel path drifted from per-node path by {difference}" diff --git a/tests/test_numba_kernel.py b/tests/test_numba_kernel.py new file mode 100644 index 0000000..812abaf --- /dev/null +++ b/tests/test_numba_kernel.py @@ -0,0 +1,179 @@ +"""The optional numba kernel must select exactly the nodes the NumPy path selects. + +``pip install numba`` swaps a compiled kernel into the best-matching-unit search. It is a second +implementation of the hottest code in the package, and the whole reason that is acceptable is that +the NumPy path stays the reference and this file asserts the two agree. + +Agreement here is **identical indices**, not close ones. Both compute +``||w||^2 - 2 x.w`` over the same centred arrays, so there is no reason for them to differ, and a +tolerance would hide the case where one of them is wrong. + +Skipped wholesale without numba. The CI job that installs it is what stops this file silently +skipping everywhere, which is the failure mode a guarded test file has. +""" + +from __future__ import annotations + +import numpy as np +import pytest + +import python_som +import python_som._som +from python_som._accelerate import bmu_kernel +from python_som._core._distance import euclidean_distance +from python_som._core._match import accumulate, bmu_indices + +BMU_KERNEL = bmu_kernel() + +pytestmark = pytest.mark.skipif(BMU_KERNEL is None, reason="needs numba") + +#: Fixed so a failure is reproducible. +SEED = 20260730 + + +def _case( + shape: tuple[int, int], n_samples: int, n_features: int, offset: float = 0.0 +) -> tuple[np.ndarray, np.ndarray]: + """Build models and a dataset, optionally far from the origin. + + :param shape: Grid shape. + :param n_samples: Number of samples. + :param n_features: Number of features. + :param offset: Constant added to both. + :return: Models and dataset. + """ + rng = np.random.default_rng(SEED) + return ( + rng.normal(size=(*shape, n_features)) + offset, + rng.normal(size=(n_samples, n_features)) + offset, + ) + + +@pytest.mark.parametrize( + ("shape", "n_samples", "n_features"), + [((5, 5), 40, 3), ((20, 20), 200, 4), ((40, 30), 500, 8), ((1, 9), 30, 2), ((60, 60), 300, 12)], +) +def test_the_kernel_selects_the_same_nodes_as_numpy( + shape: tuple[int, int], n_samples: int, n_features: int +) -> None: + """The claim the extra rests on.""" + weights, data = _case(shape, n_samples, n_features) + np.testing.assert_array_equal( + bmu_indices(data, weights, euclidean_distance, BMU_KERNEL), + bmu_indices(data, weights, euclidean_distance), + ) + + +@pytest.mark.parametrize("offset", [0.0, 1e6, 1e9, 1e12]) +def test_the_kernel_is_exact_far_from_the_origin(offset: float) -> None: + """The centring happens before the kernel is called, so it must inherit the fix. + + The kernel receives arrays that are already shifted. If it were ever changed to take raw models + and centre them itself, or not to centre at all, this is what would catch it. + """ + weights, data = _case((40, 40), 300, 6, offset=offset) + flat = weights.reshape(-1, 6) + exact = np.array([np.linalg.norm(x - flat, axis=-1).argmin() for x in data]) + np.testing.assert_array_equal(bmu_indices(data, weights, euclidean_distance, BMU_KERNEL), exact) + + +def test_the_kernel_breaks_ties_to_the_lowest_index() -> None: + """``argmin`` keeps the first minimum, and a ``<`` comparison in the kernel must match it. + + A ``<=`` in the inner loop would keep the *last* tied node instead, which no other test here + would notice. + """ + weights = np.full((3, 3, 2), 10.0) + weights[0, 2] = [1.0, 0.0] + weights[2, 0] = [1.0, 0.0] + data = np.array([[1.0, 0.0]]) + assert int(bmu_indices(data, weights, euclidean_distance, BMU_KERNEL)[0]) == 2 + + +def test_accumulate_agrees_through_the_kernel() -> None: + """Eq. (8)'s inputs must not depend on which search produced the nodes.""" + shape = (12, 9) + weights, data = _case(shape, 200, 5) + + fast_sums, fast_counts = accumulate(data, weights, shape, euclidean_distance, BMU_KERNEL) + slow_sums, slow_counts = accumulate(data, weights, shape, euclidean_distance) + + np.testing.assert_array_equal(fast_counts, slow_counts) + np.testing.assert_array_equal(fast_sums, slow_sums) + + +def test_training_agrees_end_to_end(monkeypatch: pytest.MonkeyPatch) -> None: + """A whole run, so any per-iteration divergence accumulates into view. + + Bit-identical: the two searches pick the same nodes, so Eq. (8) receives the same inputs and the + arithmetic after that point is the same code. + + The NumPy arm is produced by patching the resolver rather than by uninstalling the extra. + """ + shape, n_iteration = (20, 16), 25 + rng = np.random.default_rng(SEED) + data = rng.normal(size=(300, 6)) + initial = rng.normal(size=(*shape, 6)) + + def train() -> np.ndarray: + """Train one map with whichever backend ``_som.BMU_KERNEL`` currently names. + + :return: The trained models. + """ + som = python_som.SOM( + x=shape[0], y=shape[1], input_len=6, neighborhood_radius=3.0, random_seed=SEED + ) + som._weights = initial.copy() + som.train(data, n_iteration=n_iteration, mode="batch") + weights: np.ndarray = som.get_weights() + return weights + + accelerated = train() + monkeypatch.setattr(python_som._som, "bmu_kernel", lambda: None) + np.testing.assert_array_equal(accelerated, train()) + + +def test_a_custom_distance_still_bypasses_the_kernel() -> None: + """The kernel computes a Euclidean criterion, so it must not be reached for anything else.""" + + def manhattan(x: object, weights: object) -> np.ndarray: + """Sum of absolute differences along the last axis. + + :param x: Input vector. + :param weights: One model or an array of them. + :return: Distances. + """ + result: np.ndarray = np.abs(np.asarray(x) - np.asarray(weights)).sum(axis=-1) + return result + + weights = np.array([[[0.9, 0.9], [0.0, 1.4]]]) + data = np.array([[0.0, 0.0]]) + assert int(bmu_indices(data, weights, manhattan, BMU_KERNEL)[0]) == 1 + + +def test_importing_the_package_does_not_import_numba() -> None: + """Installing numba must not add 104 ms to every ``import python_som``. + + numba is deferred to the first call that needs it, where the JIT compile is paid anyway. A + module-level import in ``_accelerate`` would be invisible in every other test here and would + quietly undo part of what numba buys. + + A subprocess, because numba is certainly already imported in this one. + """ + import subprocess # noqa: PLC0415 + import sys # noqa: PLC0415 + + code = ( + "import sys, python_som; " + "assert 'numba' not in sys.modules, sorted(m for m in sys.modules if 'numba' in m)[:3]; " + "import numpy as np; " + "som = python_som.SOM(x=4, y=4, input_len=2, random_seed=0); " + "som.train(np.zeros((5, 2)), n_iteration=1, mode='batch'); " + "assert 'numba' in sys.modules, 'the kernel should have loaded by now'; " + "print('clean')" + ) + result = subprocess.run( # noqa: S603 + [sys.executable, "-c", code], capture_output=True, text=True, check=False + ) + assert result.returncode == 0, result.stdout + result.stderr + assert "clean" in result.stdout diff --git a/uv.lock b/uv.lock index cf37b41..1f9be63 100644 --- a/uv.lock +++ b/uv.lock @@ -2650,7 +2650,7 @@ wheels = [ [[package]] name = "python-som" -version = "0.6.1" +version = "0.7.0" source = { editable = "." } dependencies = [ { name = "numpy", version = "2.2.6", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.11'" },