diff --git a/CHANGELOG.md b/CHANGELOG.md
index e2c12b0..11c5f83 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -14,7 +14,7 @@ adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
   and
   the last training report. No pickle on either side, so loading one cannot execute code.
 
-  A reloaded map **continues the same random stream**, so training it further gives what never
+  A reloaded map continues the same random stream, so training it further gives what never
   stopping would have given. Both the seed and the generator's current state are saved: the seed
   alone
   would restart the stream and a resumed run would quietly diverge from an uninterrupted one.
@@ -31,13 +31,13 @@ adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
   and
   the member list, JSON, format version and weight shape are all validated before a map is built.
   What
-  that buys is "cannot execute code", not "safe to load anything" — an `.npz` is a zip, so a hostile
+  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. `docs/save-and-load.md` states both halves.
 
 - **`SOM.last_report`**, a frozen `TrainingReport`: mode, iterations, sample count, seed, the final
   learning rate and radius, quantization error, the `python_som` and NumPy versions, and wall time.
   `train()` still returns the quantization error, so nothing that used the return value changes.
-  `final_learning_rate` is `None` for batch training, which has no step size — Eq. (8) is a weighted
+  `final_learning_rate` is `None` for batch training, which has no step size: Eq. (8) is a weighted
   mean, and reporting the unused initial value would read as though it had been applied.
 
 - **`SOMConfig`**, the serialisable description of a map, available as `som.config()`.
@@ -85,7 +85,7 @@ adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 - The package is now a pure functional core with a thin shell around it. `python_som._core` holds
   every numeric decision as functions over NumPy arrays and imports nothing but NumPy;
   `python_som._convert` adapts pandas at the boundary; `python_som._som` keeps the state, the
-  validation and the training loops. **No public name, signature or numerical result changes** —
+  validation and the training loops. No public name, signature or numerical result changes:
   trained weights are bit-identical to 0.3.0 for every combination of training mode and neighborhood
   function, which is asserted rather than assumed.
 
@@ -94,7 +94,7 @@ adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
   point of violation. `architecture-profile.toml` records the style.
 
 - The two update rules are now pure functions returning new arrays. An earlier note claimed this
-  was also faster; **it is not**, and the claim is withdrawn. Measured with interleaved arms and
+  was also faster. It is not, and the claim is withdrawn. Measured with interleaved arms and
   medians, the pure form is about 10% slower on small maps and indistinguishable above roughly
   100x100. The cost buys functions testable without constructing a network, and is single-digit
   milliseconds over a 10,000-iteration run on a 20x20 map. `benchmarks/bench_update.py` is the
@@ -104,19 +104,19 @@ adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 - **Batch training is 1.2x to 1.5x faster**, with identical results. Eq. (8) needs the neighborhood
   between every pair of nodes, and the previous implementation got it by evaluating the neighborhood
-  function once per node, once per iteration — 1600 full-grid evaluations per iteration on a 40x40
-  map. Profiling made that **42% of batch training**, more than the contraction it feeds.
+  function once per node, once per iteration, which is 1600 full-grid evaluations per iteration on a 40x40
+  map. Profiling made that 42% of batch training, more than the contraction it feeds.
 
   It is unnecessary, because a neighborhood depends only on the offset between two nodes and never
   on where the winner sits. One kernel over every reachable offset serves the whole grid, and
-  each node's neighborhood is a slice of it — a view, not a copy. The kernel costs `4xy` floats,
+  each node's neighborhood is a slice of it, as a view rather than a copy. The kernel costs `4xy` floats,
   111 KB on a 60x60 map, against the 800 MB a full `(x, y, x, y)` tensor would need.
 
   The offset-only dependence holds on a torus as well as a flat grid, because making the fold depend
   on the offset alone is exactly what the minimum-image convention does. It holds per axis, so mixed
   maps (one axis wrapping, one not) need no special case.
 
-  **Results do not change.** Trained weights are bit-identical for all 72 working combinations of
+  Results do not change: trained weights are bit-identical for all 72 working combinations of
   training mode, neighborhood function, initializer and cyclic setting, and the kernel is checked
   against per-node evaluation across 40,832 combinations at exactly `0.0`. The gaussian gains more
   than the bubble, whose per-node form is cheaper to begin with. `benchmarks/bench_batch.py` is the
@@ -129,16 +129,16 @@ adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 - **Linear initialization was inaccurate for data far from the origin.** Through 0.3.0 its PCA went
   to scikit-learn's `auto` solver, which since 1.5 selects `covariance_eigh` when samples outnumber
-  features — it forms the covariance matrix, and squaring the data squares its condition number. On
-  `(150, 4)` data offset by 1e7, the second explained variance was wrong by **5.8%**, and the
+  features. It forms the covariance matrix, and squaring the data squares its condition number. On
+  `(150, 4)` data offset by 1e7, the second explained variance was wrong by 5.8%, and the
   resulting models differed from the correct ones by 2.43 against a total model spread of 2.0: an
   error larger than the structure being initialized.
 
   The NumPy implementation decomposes the centred matrix directly and agrees with a `longdouble`
-  reference to ~1e-15. Data offset far from the origin is not exotic — timestamps, easting/northing
+  reference to ~1e-15. Data offset far from the origin is not exotic: timestamps, easting/northing
   coordinates and absolute sensor readings all look like it.
 
-  **This changes results.** For data near the origin the difference is floating-point noise (~1e-15,
+  This changes results. For data near the origin the difference is floating-point noise (~1e-15,
   and `auto_dimensions` is unaffected because it standardizes first). Far from the origin it is
   large, and 0.4.0 is the correct one.
 
@@ -148,11 +148,11 @@ adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 ### Removed
 
 - **pandas and scikit-learn are no longer required.** `numpy` is the only runtime dependency.
-  Measured on Linux/CPython 3.12, a fresh install drops from **333 MB across 10 packages to 69 MB
-  across 1** — a 264 MB reduction, 79% of the payload, since the two also pulled in scipy, joblib,
+  Measured on Linux/CPython 3.12, a fresh install drops from 333 MB across 10 packages to 69 MB
+  across 1, a 264 MB reduction and 79% of the payload, since the two also pulled in scipy, joblib,
   narwhals, python-dateutil, six and threadpoolctl.
 
-  Nothing is lost, and input support is *wider*. pandas was used for one `isinstance` check before
+  Nothing is lost, and input support is wider. pandas was used for one `isinstance` check before
   calling `.to_numpy()`; `np.asarray` already does that via the `__array__` protocol. Because
   that is a protocol rather than a library, polars, pyarrow, xarray and CuPy objects now work too,
   without this package knowing they exist.
diff --git a/README.md b/README.md
index 262545e..7216b65 100644
--- a/README.md
+++ b/README.md
@@ -5,8 +5,8 @@
 [![Python versions](https://img.shields.io/pypi/pyversions/python-som.svg)](https://pypi.org/project/python-som/)
 [![License: MIT](https://img.shields.io/badge/License-MIT-blue.svg)](https://github.com/andremsouza/python-som/blob/master/LICENSE)
 
-Implementation of Kohonen's 2-D self-organizing map, built on NumPy, with scikit-learn for PCA and
-standardization. Accepts NumPy arrays, pandas DataFrames and plain lists.
+Implementation of Kohonen's 2-D self-organizing map. NumPy is the only dependency. Accepts NumPy
+arrays, pandas DataFrames, polars, pyarrow, and anything else implementing the `__array__` protocol.
 
 [Documentation](https://andremsouza.github.io/python-som/) ·
 [Changelog](https://github.com/andremsouza/python-som/blob/master/CHANGELOG.md)
@@ -36,7 +36,7 @@ winner = som.winner(data[0])
 ```
 
 A full worked example with plots is in [examples/iris.py](https://github.com/andremsouza/python-som/blob/master/examples/iris.py) and in the
-[getting-started guide](https://andremsouza.github.io/python-som/getting-started/).
+[getting-started guide](https://andremsouza.github.io/python-som/tutorial/first-map/).
 
 ![U-matrix of a SOM trained on Iris](https://raw.githubusercontent.com/andremsouza/python-som/master/docs/assets/iris.png)
 
@@ -68,7 +68,7 @@ update of Kohonen Eq. (8) is a weighted mean whose denominator is not sign-defin
 neighborhood function. Use `mode='random'` or `mode='sequential'`; `mode='batch'` raises a
 `ValueError`.
 
-See [Neighborhood functions](https://andremsouza.github.io/python-som/neighborhood-functions/) for
+See [Neighborhood functions](https://andremsouza.github.io/python-som/reference/neighborhood-functions/) for
 the derivations, including why the Mexican hat is not an outer product of two 1-D wavelets.
 
 ## Upgrading
diff --git a/architecture-profile.toml b/architecture-profile.toml
index b309097..60ad2bd 100644
--- a/architecture-profile.toml
+++ b/architecture-profile.toml
@@ -4,7 +4,7 @@
 # around it (`_som.py`, `_convert.py`) holds state, validation, progress output, and the pandas
 # adapter. See docs/ and `src/python_som/_core/__init__.py` for the rule in prose.
 #
-# The rule that matters day to day — the core imports nothing but NumPy — is enforced in CI by
+# The rule that matters day to day (the core imports nothing but NumPy) is enforced in CI by
 # ruff's TID251 (see `[tool.ruff.lint.flake8-tidy-imports.banned-api]` in pyproject.toml) and by
 # tests/test_core_boundary.py. This file is not run in CI: the checker that reads it is a
 # proprietary local review aid, not redistributable, so it cannot be a required check on a public
diff --git a/docs/explanation/artifact-safety.md b/docs/explanation/artifact-safety.md
new file mode 100644
index 0000000..c413dae
--- /dev/null
+++ b/docs/explanation/artifact-safety.md
@@ -0,0 +1,54 @@
+# Artifact safety
+
+Loading a `.npz` written by `save_npz` cannot execute code. Loading a `.pkl` can. That is a
+difference in kind rather than degree, and it is the reason this package ships an artifact format
+at all: not to forbid `pickle`, but to make the safe path the obvious one so nobody has to reach
+for the unsafe one.
+
+What follows is the precise boundary, because "safer" without a boundary is not a claim anyone can
+act on.
+
+Safer than a pickle, and the difference is categorical rather than a matter of degree, but not
+unconditionally safe, so here is exactly what is and is not guaranteed.
+
+**Cannot happen:**
+
+- **Code execution.** `allow_pickle=False` is passed explicitly, so a crafted object array is
+  refused by NumPy, not unpickled. The test suite builds a payload that really would create a file,
+  proves it executes when NumPy is *allowed* to unpickle it, then shows the loader refusing the same
+  file with nothing created.
+- **Loading an arbitrary function by name.** Names resolve only through the registries. Nothing from
+  the file is imported, evaluated, or looked up in a module.
+- **A partly-read file being treated as valid.** The member list, the metadata JSON, the format
+  version, the weight dimensions and the weight shape against the saved config are all checked
+  before a map is constructed. Every failure raises `ArtifactError`.
+
+**Can still happen:**
+
+- **Resource exhaustion.** An `.npz` is a zip, so a hostile file can be built to decompress to
+  something enormous. Nothing here caps that.
+
+So: treat an artifact from an untrusted source the way you would a JPEG from one. It cannot run
+code, but it can still be malformed or oversized. Do not treat it as you would a signed archive.
+
+If you need integrity as well, hash the file and record the digest alongside whatever references it;
+that is outside what this format tries to do.
+
+## Why the guarantee stops there
+
+`allow_pickle=False` closes the code-execution path and nothing else. A zip archive can still be
+crafted to decompress to something enormous, and nothing in this format caps that. Capping it
+would mean inventing a size limit with no principled value, which is the kind of unsourced
+constant this package removes when it finds one.
+
+If you need integrity as well as safety, that is a different property and wants a different tool:
+hash the file and record the digest wherever the artifact is referenced.
+
+## How the claim is tested
+
+The test suite does not assert that pickles are refused. It builds a payload with a `__reduce__`
+that creates a file, proves the payload really executes when something unpickles it, and then
+shows the loader refusing a file containing it with nothing created.
+
+Proving liveness first is the part that matters. Without it the test would pass just as well
+against an inert payload, and would be evidence of nothing.
diff --git a/docs/training-modes.md b/docs/explanation/batch-vs-stepwise.md
similarity index 97%
rename from docs/training-modes.md
rename to docs/explanation/batch-vs-stepwise.md
index a7df95a..f09ced7 100644
--- a/docs/training-modes.md
+++ b/docs/explanation/batch-vs-stepwise.md
@@ -76,7 +76,7 @@ stepwise modes. Kohonen Section 4.1 on the radius schedule:
 
 So a sensible starting point is `neighborhood_radius` at roughly half the shorter side of the grid,
 decaying from there, with the floor described in
-[Neighborhood functions](neighborhood-functions.md#the-radius-floor) keeping it away from zero.
+[Neighborhood functions](../reference/neighborhood-functions.md#the-radius-floor) keeping it away from zero.
 
 ## Reproducibility
 
diff --git a/docs/explanation/why-isotropy-matters.md b/docs/explanation/why-isotropy-matters.md
new file mode 100644
index 0000000..9ebc9d5
--- /dev/null
+++ b/docs/explanation/why-isotropy-matters.md
@@ -0,0 +1,87 @@
+# Why isotropy matters
+
+A neighborhood function has to depend on *one* number: the distance between two nodes on the grid.
+Not on the horizontal offset and the vertical offset separately. This page is about why that
+distinction is load-bearing rather than pedantic. Getting it wrong produces a function that looks
+plausible and inverts the behaviour the mexican hat exists to provide.
+
+## What the sources require
+
+Kohonen (2013) Eq. (5) writes the neighborhood as a function of `sqdist(c, i)`, "the square of the
+geometric distance between the nodes c and i in the grid". One scalar.
+
+Vrieze (1995) is more explicit still. Figure 3, *"The 'Mexican-hat' function of lateral
+interaction"*, plots interaction against an abscissa labelled **"Lateral distance"**. The text:
+
+> When a neuron is firing, near by situated neurons are stimulated, diminishing with increasing
+> distance to the firing neuron. From a certain distance on inhibition will take place that
+> gradually vanishes when the distance becomes very large.
+
+The coefficient is written $h_{i i_c} = 1/\lVert i_c - i \rVert$, the neighborhood as
+$N_i = \{i' \mid d(i,i') \le \rho\}$, and Vrieze notes that both spaces are assumed to be metric.
+Every formulation in both papers reduces the two axis offsets to a single distance *first*, then
+applies the profile.
+
+## The mistake that looks right
+
+The gaussian is **multiplicatively separable**:
+
+$$e^{-(dx^2 + dy^2)/2\sigma^2} = e^{-dx^2/2\sigma^2} \cdot e^{-dy^2/2\sigma^2}$$
+
+So you can compute it as an outer product of two 1-D gaussians and get exactly the isotropic 2-D
+gaussian. That works, it is a normal optimization, and it tempts you into thinking the same trick
+generalizes.
+
+It does not. Separability is a property of the exponential, not of neighborhood functions. Apply the
+same construction to the Ricker wavelet and you get something else entirely:
+
+```python
+ax = (1 - dx**2 / sigma**2) * exp(-(dx**2) / (2 * sigma**2))
+ay = (1 - dy**2 / sigma**2) * exp(-(dy**2) / (2 * sigma**2))
+h = outer(ax, ay)  # not a mexican hat
+```
+
+In the diagonal quadrants both 1-D factors are negative, so their product is **positive**. The
+function stops inhibiting exactly where it is supposed to inhibit most.
+
+Measured on a 21×21 grid with $\sigma = 3$, centre at (10, 10):
+
+| quantity | outer product | correct isotropic form |
+| --- | --- | --- |
+| $h$ at $(c + 2\sigma,\, c + 2\sigma)$ | +0.165 | −0.055 |
+| $h$ at $(c + 3\sigma,\, c + 3\sigma)$ | +0.008 | −0.001 |
+| global minimum | −0.443 | −0.135 (= $-e^{-2}$, at $r = 2\sigma$) |
+| zero-crossing locus | a **cross** ($dx = \pm\sigma \cup dy = \pm\sigma$) | a **circle**, $r = \sqrt{2}\sigma$ |
+
+The separable version puts an excitatory lobe worth 16.5% of the winner's strength where the
+function should be pushing models away, and along the diagonals it never becomes negative at all.
+Its value is not a function of any metric on the grid, which is precisely the property Eq. (5)
+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.
+
+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
+isotropic one.
+
+## The one deliberate exception
+
+The **bubble** is a Chebyshev ball, so it is *not* isotropic under the Euclidean metric, and that is
+intentional. Vrieze's own appendix computes `b = MAX(ABS(i - w_i), ABS(j - w_j))`, a square region
+rather than a disc. Kohonen's phrasing ("up to a certain radius from the winner") reads as Euclidean,
+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.
+
+## Further reading
+
+- [Neighborhood functions](../reference/neighborhood-functions.md) for the formulas and constants.
+- [Batch vs stepwise](batch-vs-stepwise.md) for why a signed neighborhood cannot be used with batch
+  training at all.
diff --git a/docs/how-to/reproduce-a-result.md b/docs/how-to/reproduce-a-result.md
new file mode 100644
index 0000000..73438d7
--- /dev/null
+++ b/docs/how-to/reproduce-a-result.md
@@ -0,0 +1,86 @@
+# Reproduce a result
+
+To get the same map twice, or to let someone else get the map you got.
+
+## Pin the seed
+
+```python
+som = python_som.SOM(x=10, y=10, input_len=4, random_seed=42)
+```
+
+Without `random_seed` the constructor draws one from the OS, so every run differs. The seed it used
+either way is available afterwards:
+
+```python
+som.get_random_seed()
+```
+
+The generator is per-instance. Constructing a map does not disturb NumPy's global random state, so
+two maps built in the same program do not interfere with each other and nothing else in your program
+is affected by building one.
+
+## Pin the version
+
+Numerical results are allowed to change between minor versions before 1.0.0, and between major
+versions after it. A seed alone does not pin a result across an upgrade.
+
+```
+python-som==0.4.0
+```
+
+Two specific breaks worth knowing about, 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.
+
+## Record what you ran
+
+`som.last_report` describes the run that just finished:
+
+```python
+error = som.train(data, n_iteration=100, mode=TrainingMode.BATCH)
+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',
+               numpy_version='2.5.1', wall_time_seconds=0.42)
+```
+
+That is enough to reconstruct the run: the seed, the mode, the iteration count, and the versions of
+this package and NumPy that computed it. `wall_time_seconds` is excluded from equality, so two
+identical runs produce equal reports.
+
+## Save the map itself
+
+Better than reconstructing a run is not having to:
+
+```python
+som.save_npz("figure-3.npz")
+```
+
+The file carries the models, the configuration, the seed, the generator state and the report. A
+reloaded map also continues the same random stream, so training it further gives what never stopping
+would have given.
+
+See [Save and load a map](save-and-load-a-map.md) for the details, including what happens when the
+map used a function this package cannot look up by name.
+
+## Check that you succeeded
+
+Reproducibility is a claim worth testing rather than assuming:
+
+```python
+first = build_and_train(seed=42)
+second = build_and_train(seed=42)
+assert np.abs(first.get_weights() - second.get_weights()).max() == 0.0
+```
+
+Exact equality, not `allclose`. Two runs with the same seed on the same version perform the same
+arithmetic in the same order, so any difference at all means something is not pinned that you think
+is: an unpinned dependency, a different BLAS, or data that is not identical.
diff --git a/docs/save-and-load.md b/docs/how-to/save-and-load-a-map.md
similarity index 74%
rename from docs/save-and-load.md
rename to docs/how-to/save-and-load-a-map.md
index ec49a8d..18d645d 100644
--- a/docs/save-and-load.md
+++ b/docs/how-to/save-and-load-a-map.md
@@ -74,7 +74,7 @@ python_som.SOM.load_npz("checkpoint.npz").train(data, n_iteration=100)
 ## The training report
 
 `som.last_report` describes the most recent `train()` call, and `None` before the first one. It is
-saved with the map, so a figure can be traced back to the run that produced it — the seed, the
+saved with the map, so a figure can be traced back to the run that produced it: the seed, the
 iteration count, the error, and the versions of this package and NumPy that computed it.
 
 `train()` still returns the quantization error, so nothing that used the return value changes.
@@ -100,39 +100,17 @@ python_som.SOM.load_npz("custom.npz", distance_function=cosine)  # works
 ```
 
 The same applies to a `functools.partial`: it has no name to record, so it is treated as custom.
-That is the correct outcome — resolving `partial(exponential_decay, factor=3.0)` by the wrapped
+That is the correct outcome. Resolving `partial(exponential_decay, factor=3.0)` by the wrapped
 function's name would silently restore a different decay from the one that trained the map.
 
 The names that *do* resolve are the keys of `NEIGHBORHOOD_FUNCTIONS`, `DECAY_FUNCTIONS` and
 `DISTANCE_FUNCTIONS`, each of which is simply the function's own name.
 
-## How safe is it to load a file from someone else
+## Is it safe to load a file from someone else
 
-Safer than a pickle, and the difference is categorical rather than a matter of degree — but not
-unconditionally safe, so here is exactly what is and is not guaranteed.
-
-**Cannot happen:**
-
-- **Code execution.** `allow_pickle=False` is passed explicitly, so a crafted object array is
-  refused by NumPy, not unpickled. The test suite builds a payload that really would create a file,
-  proves it executes when NumPy is *allowed* to unpickle it, then shows the loader refusing the same
-  file with nothing created.
-- **Loading an arbitrary function by name.** Names resolve only through the registries. Nothing from
-  the file is imported, evaluated, or looked up in a module.
-- **A partly-read file being treated as valid.** The member list, the metadata JSON, the format
-  version, the weight dimensions and the weight shape against the saved config are all checked
-  before a map is constructed. Every failure raises `ArtifactError`.
-
-**Can still happen:**
-
-- **Resource exhaustion.** An `.npz` is a zip, so a hostile file can be built to decompress to
-  something enormous. Nothing here caps that.
-
-So: treat an artifact from an untrusted source the way you would a JPEG from one — it cannot run
-code, but it can still be malformed or oversized. Do not treat it as you would a signed archive.
-
-If you need integrity as well, hash the file and record the digest alongside whatever references it;
-that is outside what this format tries to do.
+Safer than a pickle, and categorically so: an artifact cannot execute code. It can still be
+malformed or oversized. [Artifact safety](../explanation/artifact-safety.md) sets out exactly what is
+and is not guaranteed, and why the guarantee stops where it does.
 
 ## Version differences
 
diff --git a/docs/how-to/use-a-custom-strategy.md b/docs/how-to/use-a-custom-strategy.md
new file mode 100644
index 0000000..f2f21ce
--- /dev/null
+++ b/docs/how-to/use-a-custom-strategy.md
@@ -0,0 +1,56 @@
+# Use a custom strategy
+
+Three things can be replaced with your own implementation: the neighborhood, the decay applied to
+the learning rate and radius, and the distance.
+
+Three things can be replaced with your own implementation: the neighborhood, the decay applied to
+the learning rate and radius, and the distance. Each has a `Protocol` describing what it must accept
+and return, so a type checker verifies yours against a named contract rather than a bare `Callable`.
+
+```python
+from python_som import DecayFunction
+
+
+def linear_to_zero(value: float, step: int, total: int) -> float:
+    return value * (1.0 - step / total)
+
+
+som = python_som.SOM(x=10, y=10, input_len=4, learning_rate_decay=linear_to_zero)
+```
+
+The protocols are **structural**: nothing needs to inherit from them, and every callable that
+already worked still does. Their parameters are positional-only, so your parameter names are your
+own. `linear_to_zero(value, step, total)` and `linear_to_zero(a, b, c)` both satisfy
+`DecayFunction`.
+
+One contract is worth reading before writing a neighborhood of your own:
+[`NeighborhoodFunction`][python_som.NeighborhoodFunction] must be a function of the grid distance
+between two nodes, not of the two axis offsets separately. See
+[Why isotropy matters](../explanation/why-isotropy-matters.md) for what goes wrong otherwise.
+
+## 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.
+
+```python
+python_som.SOM.load_npz("custom.npz")
+# ArtifactError: ... Pass them back explicitly: load_npz(path, distance_function=...)
+
+python_som.SOM.load_npz("custom.npz", distance_function=my_distance)  # works
+```
+
+That is deliberate. The alternative is guessing, and a map that silently reloads with a different
+distance from the one that trained it is worse than one that refuses.
+
+If you want a strategy that round-trips without the extra argument, the ones this package ships all
+do: their names are the keys of `NEIGHBORHOOD_FUNCTIONS`, `DECAY_FUNCTIONS` and
+`DISTANCE_FUNCTIONS`.
+
+## Before writing a neighborhood
+
+Read [Why isotropy matters](../explanation/why-isotropy-matters.md) first. A neighborhood must be a
+function of the grid distance between two nodes, not of the two axis offsets separately, and the
+construction that tempts you into getting this wrong is the same one that works correctly for the
+gaussian.
diff --git a/docs/api.md b/docs/reference/api.md
similarity index 100%
rename from docs/api.md
rename to docs/reference/api.md
diff --git a/docs/changelog.md b/docs/reference/changelog.md
similarity index 100%
rename from docs/changelog.md
rename to docs/reference/changelog.md
diff --git a/docs/neighborhood-functions.md b/docs/reference/neighborhood-functions.md
similarity index 100%
rename from docs/neighborhood-functions.md
rename to docs/reference/neighborhood-functions.md
diff --git a/docs/options-and-types.md b/docs/reference/options.md
similarity index 71%
rename from docs/options-and-types.md
rename to docs/reference/options.md
index fbbe5f9..1c5f305 100644
--- a/docs/options-and-types.md
+++ b/docs/reference/options.md
@@ -26,7 +26,7 @@ before stops accepting one.
 ## Why bother
 
 A type checker cannot tell `"batch"` from `"bacth"`, so the misspelling survives until the
-`ValueError` at runtime — after however long it took to get there. The parameters are typed as the
+`ValueError` at runtime, after however long it took to get there. The parameters are typed as the
 enum *or* the exact set of valid strings, so both spellings pass and a typo does not:
 
 ```python
@@ -73,39 +73,12 @@ Strings are deprecated as of **0.4.0**, and nothing warns yet.
 No warning in 0.4.0 is deliberate. `mode="batch"` is what this documentation showed until this page
 existed, so a minor release that warned on it would be scolding people for following its own
 instructions. The written notice comes first and the warning follows in 0.5.0, which also means
-0.5.0 has to ship before 1.0.0 can remove anything — the policy is that a removal in a major release
+0.5.0 has to ship before 1.0.0 can remove anything: the policy is that a removal in a major release
 is preceded by at least one minor release that warns.
 
 Migrating early costs nothing and is a mechanical substitution: `"batch"` becomes
 `TrainingMode.BATCH`, and so on down the table above.
 
-## Custom strategies
-
-Three things can be replaced with your own implementation: the neighborhood, the decay applied to
-the learning rate and radius, and the distance. Each has a `Protocol` describing what it must accept
-and return, so a type checker verifies yours against a named contract rather than a bare `Callable`.
-
-```python
-from python_som import DecayFunction
-
-
-def linear_to_zero(value: float, step: int, total: int) -> float:
-    return value * (1.0 - step / total)
-
-
-som = python_som.SOM(x=10, y=10, input_len=4, learning_rate_decay=linear_to_zero)
-```
-
-The protocols are **structural**: nothing needs to inherit from them, and every callable that
-already worked still does. Their parameters are positional-only, so your parameter names are your
-own — `linear_to_zero(value, step, total)` and `linear_to_zero(a, b, c)` both satisfy
-`DecayFunction`.
-
-One contract is worth reading before writing a neighborhood of your own:
-[`NeighborhoodFunction`][python_som.NeighborhoodFunction] must be a function of the grid distance
-between two nodes, not of the two axis offsets separately. See
-[Neighborhood functions](neighborhood-functions.md) for what goes wrong otherwise.
-
 ## Learning rate
 
 `learning_rate` is validated from 0.4.0, having been unchecked before.
@@ -114,7 +87,7 @@ between two nodes, not of the two axis offsets separately. See
   completes and changes nothing; `-1` moves models *away* from the samples they match, taking the
   quantization error from 0.0 to 11.7 in one run. Both used to be accepted in silence.
 - **Warned about**: anything above 1. Eq. (3) moves a model a fraction `alpha * h` of the way to the
-  sample, so above 1 it overshoots and oscillates. It does not necessarily diverge — at `alpha = 5`
+  sample, so above 1 it overshoots and oscillates. It does not necessarily diverge: at `alpha = 5`
   with decay disabled the largest weight stayed at 3.61, because the neighborhood damps the
-  correction away from the winner — and Kohonen gives no upper bound, so it is a warning rather than
+  correction away from the winner. Kohonen gives no upper bound, so it is a warning rather than
   an error.
diff --git a/docs/getting-started.md b/docs/tutorial/first-map.md
similarity index 93%
rename from docs/getting-started.md
rename to docs/tutorial/first-map.md
index 1b625e6..ed906a5 100644
--- a/docs/getting-started.md
+++ b/docs/tutorial/first-map.md
@@ -73,7 +73,7 @@ error = som.train(features, n_iteration=len(features), mode="batch", verbose=Tru
 print(f"quantization error: {error:.4f}")
 ```
 
-See [Training modes](training-modes.md) for how the three modes differ and why batch is
+See [Batch vs stepwise](../explanation/batch-vs-stepwise.md) for how the modes differ and why batch is
 recommended.
 
 ### Inspect the result
@@ -109,7 +109,7 @@ plt.axis((0, som.get_shape()[0], 0, som.get_shape()[1]))
 plt.savefig("iris.png")
 ```
 
-![U-matrix of a SOM trained on Iris](assets/iris.png)
+![U-matrix of a SOM trained on Iris](../assets/iris.png)
 
 Darker regions of the U-matrix mean neighbouring models are far apart, so they read as boundaries
 between clusters.
@@ -135,4 +135,4 @@ print(som.get_random_seed())
 
 Same seed, same map. The generator belongs to the instance, so building a SOM will not disturb
 NumPy's global random state. Note that seeds do not reproduce across the 0.3.0 boundary; see the
-warning in [Training modes](training-modes.md#reproducibility).
+warning in [Batch vs stepwise](../explanation/batch-vs-stepwise.md#reproducibility).
diff --git a/mkdocs.yml b/mkdocs.yml
index 4758327..1594041 100644
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -29,16 +29,35 @@ theme:
 
 nav:
   - Home: index.md
-  - Getting started: getting-started.md
-  - Neighborhood functions: neighborhood-functions.md
-  - Training modes: training-modes.md
-  - Options and types: options-and-types.md
-  - Saving and loading: save-and-load.md
-  - API reference: api.md
-  - Changelog: changelog.md
+  - Tutorial:
+      - Your first map: tutorial/first-map.md
+  - How-to:
+      - Save and load a map: how-to/save-and-load-a-map.md
+      - Reproduce a result: how-to/reproduce-a-result.md
+      - Use a custom strategy: how-to/use-a-custom-strategy.md
+  - Reference:
+      - Options and types: reference/options.md
+      - Neighborhood functions: reference/neighborhood-functions.md
+      - API: reference/api.md
+      - Changelog: reference/changelog.md
+  - Explanation:
+      - Why isotropy matters: explanation/why-isotropy-matters.md
+      - Batch vs stepwise: explanation/batch-vs-stepwise.md
+      - Artifact safety: explanation/artifact-safety.md
 
 plugins:
   - search
+  # Every page moved when the docs were reorganised, and two of the old URLs are frozen in 0.3.0's
+  # PyPI description, which can never be edited. These keep them working.
+  - redirects:
+      redirect_maps:
+        getting-started.md: tutorial/first-map.md
+        neighborhood-functions.md: reference/neighborhood-functions.md
+        training-modes.md: explanation/batch-vs-stepwise.md
+        options-and-types.md: reference/options.md
+        save-and-load.md: how-to/save-and-load-a-map.md
+        api.md: reference/api.md
+        changelog.md: reference/changelog.md
   - mkdocstrings:
       handlers:
         python:
diff --git a/pyproject.toml b/pyproject.toml
index 070f4a0..d9a2229 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -66,6 +66,16 @@ 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.
+    "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
diff --git a/uv.lock b/uv.lock
index 77c3222..61d7020 100644
--- a/uv.lock
+++ b/uv.lock
@@ -1606,6 +1606,18 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/5b/54/662a4743aa81d9582ee9339d4ffa3c8fd40a4965e033d77b9da9774d3960/mkdocs_material_extensions-1.3.1-py3-none-any.whl", hash = "sha256:adff8b62700b25cb77b53358dad940f3ef973dd6db797907c49e3c2ef3ab4e31", size = 8728, upload-time = "2023-11-22T19:09:43.465Z" },
 ]
 
+[[package]]
+name = "mkdocs-redirects"
+version = "1.2.2"
+source = { registry = "https://pypi.org/simple" }
+dependencies = [
+    { name = "mkdocs" },
+]
+sdist = { url = "https://files.pythonhosted.org/packages/f1/a8/6d44a6cf07e969c7420cb36ab287b0669da636a2044de38a7d2208d5a758/mkdocs_redirects-1.2.2.tar.gz", hash = "sha256:3094981b42ffab29313c2c1b8ac3969861109f58b2dd58c45fc81cd44bfa0095", size = 7162, upload-time = "2024-11-07T14:57:21.109Z" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/c4/ec/38443b1f2a3821bbcb24e46cd8ba979154417794d54baf949fefde1c2146/mkdocs_redirects-1.2.2-py3-none-any.whl", hash = "sha256:7dbfa5647b79a3589da4401403d69494bd1f4ad03b9c15136720367e1f340ed5", size = 6142, upload-time = "2024-11-07T14:57:19.143Z" },
+]
+
 [[package]]
 name = "mkdocstrings"
 version = "1.0.6"
@@ -2582,6 +2594,7 @@ dev = [
 ]
 docs = [
     { name = "mkdocs-material" },
+    { name = "mkdocs-redirects" },
     { name = "mkdocstrings-python" },
 ]
 examples = [
@@ -2600,6 +2613,7 @@ requires-dist = [
     { name = "hypothesis", marker = "extra == 'dev'", specifier = "==6.163.0" },
     { name = "matplotlib", marker = "extra == 'examples'", specifier = ">=3.8" },
     { name = "mkdocs-material", marker = "extra == 'docs'", specifier = "==9.7.7" },
+    { name = "mkdocs-redirects", marker = "extra == 'docs'", specifier = "==1.2.2" },
     { name = "mkdocstrings-python", marker = "extra == 'docs'", specifier = "==2.0.5" },
     { name = "mypy", marker = "extra == 'dev'", specifier = "==2.3.0" },
     { name = "numpy", specifier = ">=1.24" },