diff --git a/README.md b/README.md
index 93496b3..70183f8 100644
--- a/README.md
+++ b/README.md
@@ -11,7 +11,7 @@ Most features were implemented using NumPy, with Scikit-learn for standardizatio
 * Linear weight initialization (with PCA)
 * Automatic selection of map size ratio (with PCA)
 * Support for cyclic arrays, for toroidal or spherical maps
-* Gaussian and Bubble neighborhood functions
+* 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)
@@ -108,6 +108,22 @@ The following are lists of public methods and functions currently available in t
 * SOM.train
 * SOM.weight_initialization
 
+## Neighborhood functions
+The `neighborhood_function` parameter selects how the winner's correction is spread over the grid.
+All three are functions of the distance between two nodes in the grid, `sqdist(c, i)` in Eq. (5) of
+Kohonen (2013):
+
+| Name | Shape | Notes |
+| --- | --- | --- |
+| `'gaussian'` | `exp(-r² / 2σ²)` | Strictly positive and monotonically decreasing. The default. |
+| `'bubble'` | `1` for `r ≤ σ`, else `0` | The truncated inner lobe of the mexican hat. Vrieze (1995) notes this flat choice is "just as effective and sometimes even better". |
+| `'mexicanhat'` | `(1 - u)·exp(-u)`, `u = r² / 2σ²` | Excitatory near the winner, inhibitory beyond it. Crosses zero at `r = √2·σ`, reaches its minimum of `-e⁻² ≈ -0.135` at `r = 2σ`. |
+
+The mexican hat takes negative values, so it **cannot be used with `mode='batch'`**: the batch update
+of Kohonen (2013), Eq. (8), is a weighted mean whose denominator `Σⱼ nⱼ·hⱼᵢ` is not sign-definite for
+a signed neighborhood function. Use `mode='random'` or `mode='sequential'` instead; passing
+`mode='batch'` raises a `ValueError`.
+
 ## References
 This implemetation was based on the following paper, by Professor Teuvo Kohonen:
 
@@ -118,4 +134,15 @@ Volume 37,
 2013,
 Pages 52-65,
 ISSN 0893-6080,
-https://doi.org/10.1016/j.neunet.2012.09.018.
\ No newline at end of file
+https://doi.org/10.1016/j.neunet.2012.09.018.
+
+The mexican hat neighborhood function follows the lateral-interaction formulation in:
+
+O. J. Vrieze,
+Kohonen network,
+in: Artificial Neural Networks: An Introduction to ANN Theory and Practice,
+Lecture Notes in Computer Science, Volume 931,
+Springer, Berlin, Heidelberg,
+1995,
+Pages 83-100,
+https://doi.org/10.1007/BFb0027024.
\ No newline at end of file
diff --git a/python_som/__init__.py b/python_som/__init__.py
index b5ae0bd..1422ba8 100644
--- a/python_som/__init__.py
+++ b/python_som/__init__.py
@@ -9,7 +9,7 @@
     - Linear weight initialization (with PCA)
     - Automatic selection of map size ratio (with PCA)
     - Support for cyclic arrays, for toroidal or spherical maps
-    - Gaussian and Bubble neighborhood functions
+    - 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)
@@ -124,7 +124,7 @@ class SOM:
         - Linear weight initialization (with PCA)
         - Automatic selection of map size ratio (with PCA)
         - Support for cyclic arrays, for toroidal or spherical maps
-        - Gaussian and Bubble neighborhood functions
+        - 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)
@@ -190,7 +190,7 @@ def __init__(
             variable. May be a predefined one from this package, or a custom function, with the
             same parameters and return type. Defaults to _asymptotic_decay
         :param neighborhood_function: str: Neighborhood function name for the training process.
-            May be either 'gaussian' or 'bubble'.
+            May be either 'gaussian', 'bubble' or 'mexicanhat'.
         :param distance_function: function: Function for calculating distances/dissimilarities
             between models of the network.
             May be a predefined one from this package, or a custom function, with the same
@@ -242,10 +242,19 @@ def __init__(
         self._learning_rate_decay = learning_rate_decay
         self._neighborhood_radius = float(neighborhood_radius)
         self._neighborhood_radius_decay = neighborhood_radius_decay
-        self._neighborhood_function = {
+        neighborhood_functions = {
             "gaussian": self._gaussian,
             "bubble": self._bubble,
-        }[neighborhood_function]
+            "mexicanhat": self._mexicanhat,
+        }
+        try:
+            self._neighborhood_function = neighborhood_functions[neighborhood_function]
+        except KeyError as exc:
+            raise ValueError(
+                "Invalid value for 'neighborhood_function' parameter. Value should be in "
+                + str(list(neighborhood_functions.keys()))
+            ) from exc
+        self._neighborhood_function_name = neighborhood_function
         self._distance_function = distance_function
         self._cyclic = (bool(cyclic_x), bool(cyclic_y))
         self._neigx, self._neigy = np.arange(self._shape[0]), np.arange(self._shape[1])
@@ -484,6 +493,18 @@ def train(
         elif mode == "sequential":
             self._train_sequential(data_array, n_iteration, verbose)
         elif mode == "batch":
+            if self._neighborhood_function_name == "mexicanhat":
+                # The batch update of Kohonen (2013), Eq. (8), is a weighted mean whose denominator
+                # is sum_j n_j * h_ji. That is only well defined for a non-negative neighborhood
+                # function. The mexican hat is signed by construction, so the denominator can vanish
+                # or turn negative, which makes the update blow up or invert the sign of the
+                # correction. Use a stepwise mode instead.
+                raise ValueError(
+                    "The 'mexicanhat' neighborhood function cannot be used with the 'batch'"
+                    " training mode: the weighted mean of Kohonen (2013), Eq. (8), is undefined"
+                    " for a neighborhood function that takes negative values, as its denominator"
+                    " is not sign-definite. Use mode='random' or mode='sequential' instead."
+                )
             self._train_batch(data_array, n_iteration, verbose)
         else:
             # Invalid training mode value
@@ -641,7 +662,7 @@ def _train_batch(
     def weight_initialization(
         self,
         mode: str = "random",
-        **kwargs: np.ndarray | pd.DataFrame | list | str | int
+        **kwargs: np.ndarray | pd.DataFrame | list | str | int,
     ) -> None:
         """Function for weight initialization of the self-organizing map.
 
@@ -758,28 +779,70 @@ def _weight_initialization_sample(
         )
         self._weights = data_array[sample].reshape(self._weights.shape)
 
+    @staticmethod
+    def _axis_offsets(
+        coords: np.ndarray, center: int, length: float, cyclic: bool
+    ) -> np.ndarray:
+        """
+        Signed offsets from 'center' to each coordinate along one axis of the grid.
+
+        When the axis is cyclic, the minimum-image convention is applied, i.e., each offset is
+        folded into the interval [-length/2, length/2), so that 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.
+
+        :param coords: np.ndarray: Coordinates along the axis, i.e., np.arange(length).
+        :param center: int: Coordinate of the center (winner node) along the axis.
+        :param length: float: Number of nodes along the axis.
+        :param cyclic: bool: Whether the axis wraps around.
+        :return: np.ndarray: Signed offsets from center to each coordinate.
+        """
+        d = coords - center
+        if cyclic:
+            d = (d + length / 2) % length - length / 2
+        return d
+
+    def _sqdist(self, c: tuple[int, ...], sigma: float) -> np.ndarray:
+        """
+        Squared geometric distance from the center c to every node of the grid.
+
+        This is the sqdist(c, i) term of Kohonen (2013), Eq. (5). Every neighborhood function must
+        be a function of this scalar quantity: the neighborhood describes lateral interaction over
+        a metric space, so it depends on the distance between two nodes and not on the individual
+        axis offsets. Only the gaussian is separable into a product of per-axis terms, and that is
+        a property of the exponential rather than of neighborhood functions in general.
+
+        :param c: (int, int): Center coordinates.
+        :param sigma: float: Spread variable, validated here for all neighborhood functions.
+        :return: np.ndarray: Squared distances from c to each node, with the shape of the network.
+        :raises ValueError: If sigma is not strictly positive.
+        """
+        if not np.isfinite(sigma) or sigma <= 0:
+            raise ValueError(
+                "The neighborhood radius 'sigma' must be a finite positive number, got "
+                + str(sigma)
+            )
+        dx = self._axis_offsets(
+            self._neigx, c[0], float(self._shape[0]), self._cyclic[0]
+        )
+        dy = self._axis_offsets(
+            self._neigy, c[1], float(self._shape[1]), self._cyclic[1]
+        )
+        return np.add.outer(np.power(dx, 2), np.power(dy, 2))
+
     def _gaussian(self, c: tuple[int, ...], sigma: float) -> np.ndarray:
         """
         Gaussian neighborhood function, centered in c. Has support for cyclic arrays.
 
+        Implements h_ci = exp(-sqdist(c, i) / (2 * sigma^2)), i.e., Eq. (5) of Kohonen (2013)
+        with the learning rate factored out, so that h_cc = 1.
+
         :param c: (int, int): Center coordinates for gaussian function.
-        :param sigma: float: Spread variable for gaussian function.
+        :param sigma: float: Spread variable for gaussian function. Must be positive.
         :return: np.ndarray: Gaussian, centered in c, over all the weights of the network.
+        :raises ValueError: If sigma is not strictly positive.
         """
-        # Calculate coefficient with sigma
-        d = 2 * sigma * sigma
-        # Calculate vertical and horizontal distances
-        dx = self._neigx - c[0]
-        dy = self._neigy - c[1]
-        # If using cyclic arrays, perform fold back distance
-        if self._cyclic[0]:
-            dx[dx > self._shape[0] / 2] -= self._shape[0]
-        if self._cyclic[1]:
-            dy[dy > self._shape[1] / 2] -= self._shape[1]
-        # Calculate gaussian centered in c
-        ax = np.exp(-np.power(dx, 2) / d)
-        ay = np.exp(-np.power(dy, 2) / d)
-        return np.outer(ax, ay)
+        return np.exp(-self._sqdist(c, sigma) / (2 * sigma * sigma))
 
     def _bubble(self, c: tuple[int, ...], sigma: float) -> np.ndarray:
         """
@@ -788,10 +851,26 @@ def _bubble(self, c: tuple[int, ...], sigma: float) -> np.ndarray:
         The neighbors of c are the nodes in the region of sigma positions in the vertical
         and horizontal directions around c.
 
+        This is the truncated inner, excitatory lobe of the mexican hat lateral-interaction
+        function; cf. 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."
+
+        Unlike the other neighborhood functions, a radius of zero is admissible here: it selects
+        the winner alone, which is well defined, whereas for the gaussian and the mexican hat a
+        zero radius is a division by zero.
+
         :param c: (int, int): Center coordinates for gaussian function.
-        :param sigma: float: Spread variable for gaussian function.
+        :param sigma: float: Spread variable for gaussian function. Must be finite and
+            non-negative.
         :return: np.ndarray: Neighborhood matrix, centered in c.
+        :raises ValueError: If sigma is negative or not finite.
         """
+        if not np.isfinite(sigma) or sigma < 0:
+            raise ValueError(
+                "The neighborhood radius 'sigma' must be a finite non-negative number, got "
+                + str(sigma)
+            )
         # Convert sigma to integer
         sigma = int(np.around(sigma))
 
@@ -813,6 +892,46 @@ def _bubble(self, c: tuple[int, ...], sigma: float) -> np.ndarray:
 
         return np.outer(ax, ay).astype(int)
 
+    def _mexicanhat(self, c: tuple[int, ...], sigma: float) -> np.ndarray:
+        """
+        Mexican Hat neighborhood function, centered in c. Has support for cyclic arrays.
+
+        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 beyond a
+        certain distance are inhibited, and the inhibition vanishes as the distance grows further
+        (Vrieze, 1995, Fig. 3).
+
+        Implements the isotropic 2-D form as a function of u = sqdist(c, i) / (2 * sigma^2):
+
+            h_ci = (1 - u) * exp(-u)
+
+        normalized so that h_cc = 1. It crosses zero at a radius of sqrt(2) * sigma, reaches its
+        minimum of -exp(-2) ~= -0.135 at a radius of 2 * sigma, and decays to zero thereafter.
+
+        Note that this is *not* the outer product of two 1-D Ricker wavelets. Unlike the gaussian,
+        the Ricker wavelet is not multiplicatively separable, and such a product is not a function
+        of the distance between two nodes: it would be positive in the diagonal quadrants, where
+        both 1-D factors are negative, placing a spurious excitatory lobe exactly where the
+        function must inhibit. The neighborhood function must depend on the grid distance alone,
+        cf. sqdist(c, i) in Kohonen (2013), Eq. (5), and the "lateral distance" abscissa of
+        Vrieze (1995), Fig. 3.
+
+        The '_bubble' neighborhood function is the truncated inner (excitatory) lobe of this same
+        lateral-interaction function; cf. Vrieze (1995), p. 85.
+
+        :param c: (int, int): Center coordinates for mexican hat function.
+        :param sigma: float: Spread variable for mexican hat function. Must be positive.
+        :return: np.ndarray: Mexican Hat, centered in c, over all the weights of the network.
+        :raises ValueError: If sigma is not strictly positive.
+
+        References:
+        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.
+        """
+        u = self._sqdist(c, sigma) / (2 * sigma * sigma)
+        return (1.0 - u) * np.exp(-u)
+
     def _data_to_numpy(self, data: np.ndarray | pd.DataFrame | list) -> np.ndarray:
         """
         Converts data to numpy array.
diff --git a/tests/test_neighborhood.py b/tests/test_neighborhood.py
new file mode 100644
index 0000000..a504c49
--- /dev/null
+++ b/tests/test_neighborhood.py
@@ -0,0 +1,290 @@
+"""Analytic property tests for the neighborhood functions of the self-organizing map.
+
+These assert closed-form properties rather than golden values, so they remain meaningful if the
+implementation is rewritten. The central property is *isotropy*: a neighborhood function models
+lateral interaction over the grid, so it must be a function of the distance between two nodes and
+nothing else.
+
+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 -- Eq. (5), the neighborhood as a function of
+sqdist(c, i).
+
+O. J. Vrieze, Kohonen network, in: Artificial Neural Networks: An Introduction to ANN Theory and
+Practice, LNCS 931, Springer, 1995, pp. 83-100, https://doi.org/10.1007/BFb0027024 -- Fig. 3, the
+"Mexican-hat" function of lateral interaction plotted against a single "lateral distance" axis.
+"""
+
+import numpy as np
+import pytest
+
+import python_som
+
+
+def make_som(
+    x=21, y=21, cyclic_x=False, cyclic_y=False, neighborhood_function="gaussian"
+):
+    """Builds a SOM with a fixed seed; only the grid geometry matters for these tests."""
+    return python_som.SOM(
+        x=x,
+        y=y,
+        input_len=3,
+        neighborhood_function=neighborhood_function,
+        cyclic_x=cyclic_x,
+        cyclic_y=cyclic_y,
+        random_seed=0,
+    )
+
+
+def grid_radius(shape, c):
+    """Euclidean distance from c to every node of a non-cyclic grid."""
+    dx = np.arange(shape[0]) - c[0]
+    dy = np.arange(shape[1]) - c[1]
+    return np.sqrt(np.add.outer(dx**2, dy**2))
+
+
+# --------------------------------------------------------------------------------------------
+# Properties shared by every neighborhood function
+# --------------------------------------------------------------------------------------------
+
+
+@pytest.mark.parametrize("name", ["gaussian", "bubble", "mexicanhat"])
+def test_shape_matches_the_grid(name):
+    som = make_som(neighborhood_function=name)
+    assert som._neighborhood_function((10, 10), 3.0).shape == (21, 21)
+
+
+@pytest.mark.parametrize("name", ["gaussian", "bubble", "mexicanhat"])
+def test_winner_has_maximum_response(name):
+    """The winner is maximally excited: h(c, c) is the global maximum."""
+    som = make_som(neighborhood_function=name)
+    h = som._neighborhood_function((10, 10), 3.0)
+    assert h[10, 10] == pytest.approx(h.max())
+
+
+@pytest.mark.parametrize("name", ["gaussian", "mexicanhat"])
+def test_normalized_to_unity_at_the_winner(name):
+    som = make_som(neighborhood_function=name)
+    assert som._neighborhood_function((10, 10), 3.0)[10, 10] == pytest.approx(1.0)
+
+
+@pytest.mark.parametrize("name", ["gaussian", "bubble", "mexicanhat"])
+@pytest.mark.parametrize("sigma", [-1.0, -0.5, float("nan"), float("inf")])
+def test_negative_or_non_finite_sigma_is_rejected(name, sigma):
+    """A negative or non-finite radius is meaningless for any neighborhood function."""
+    som = make_som(neighborhood_function=name)
+    with pytest.raises(ValueError, match="sigma"):
+        som._neighborhood_function((10, 10), sigma)
+
+
+@pytest.mark.parametrize("name", ["gaussian", "mexicanhat"])
+def test_zero_sigma_is_rejected_where_it_divides(name):
+    """sigma appears in a denominator, so zero must raise rather than yield inf/nan."""
+    som = make_som(neighborhood_function=name)
+    with pytest.raises(ValueError, match="sigma"):
+        som._neighborhood_function((10, 10), 0.0)
+
+
+def test_bubble_admits_a_zero_radius():
+    """For the bubble a zero radius is well defined: it selects the winner alone."""
+    som = make_som(neighborhood_function="bubble")
+    h = som._bubble((10, 10), 0.0)
+    assert h[10, 10] == 1
+    assert h.sum() == 1
+
+
+@pytest.mark.parametrize("name", ["gaussian", "bubble", "mexicanhat"])
+def test_four_fold_symmetry_about_the_winner(name):
+    """h is symmetric under reflection in both axes when the winner is centred."""
+    som = make_som(neighborhood_function=name)
+    h = som._neighborhood_function((10, 10), 3.0)
+    np.testing.assert_allclose(h, h[::-1, :])
+    np.testing.assert_allclose(h, h[:, ::-1])
+
+
+@pytest.mark.parametrize("name", ["gaussian", "bubble", "mexicanhat"])
+def test_is_isotropic(name):
+    """h depends only on the grid distance -- the property a separable product violates.
+
+    Kohonen (2013) Eq. (5) defines the neighborhood over sqdist(c, i); Vrieze (1995) Fig. 3 plots
+    it against a single "lateral distance" axis. Nodes equidistant from the winner must therefore
+    receive equal responses. An outer product of two 1-D profiles fails this for every profile
+    except the gaussian, which is separable only by accident of the exponential.
+    """
+    som = make_som(neighborhood_function=name)
+    h = som._neighborhood_function((10, 10), 3.0)
+    r = np.round(grid_radius((21, 21), (10, 10)), 9)
+    for radius in np.unique(r):
+        responses = h[r == radius]
+        np.testing.assert_allclose(responses, responses[0], atol=1e-12)
+
+
+# --------------------------------------------------------------------------------------------
+# Gaussian
+# --------------------------------------------------------------------------------------------
+
+
+def test_gaussian_matches_its_closed_form():
+    som = make_som()
+    sigma = 3.0
+    h = som._gaussian((10, 10), sigma)
+    expected = np.exp(-(grid_radius((21, 21), (10, 10)) ** 2) / (2 * sigma**2))
+    np.testing.assert_allclose(h, expected)
+
+
+def test_gaussian_is_strictly_positive_and_decreasing():
+    som = make_som()
+    h = som._gaussian((10, 10), 3.0)
+    assert (h > 0).all()
+    centre_row = h[10, 10:]
+    assert (np.diff(centre_row) < 0).all()
+
+
+# --------------------------------------------------------------------------------------------
+# Mexican hat -- the properties that define the shape
+# --------------------------------------------------------------------------------------------
+
+
+def test_mexican_hat_crosses_zero_at_sqrt2_sigma():
+    """(1 - u) exp(-u) with u = r^2/(2 sigma^2) vanishes at u = 1, i.e. r = sqrt(2) sigma."""
+    som = make_som(x=41, y=41)
+    sigma = 4.0
+    h = som._mexicanhat((20, 20), sigma)
+    r = grid_radius((41, 41), (20, 20))
+    inside = r < np.sqrt(2) * sigma - 1e-9
+    outside_but_near = (r > np.sqrt(2) * sigma + 1e-9) & (r < 4 * sigma)
+    assert (h[inside] > 0).all()
+    assert (h[outside_but_near] < 0).all()
+
+
+def test_mexican_hat_minimum_is_minus_exp_minus_two_at_two_sigma():
+    """d/du [(1-u) e^-u] = 0 at u = 2, giving h = -e^-2 at r = 2 sigma."""
+    som = make_som(x=41, y=41)
+    sigma = 4.0
+    h = som._mexicanhat((20, 20), sigma)
+    r = grid_radius((41, 41), (20, 20))
+    assert h.min() == pytest.approx(-np.exp(-2.0), rel=1e-9)
+    assert r.flat[h.argmin()] == pytest.approx(2 * sigma, rel=1e-9)
+
+
+def test_mexican_hat_is_inhibitory_on_the_diagonal():
+    """Regression for the separable outer product.
+
+    An outer product of two 1-D Ricker wavelets is positive wherever both factors are negative,
+    putting a spurious excitatory lobe on the diagonals. Measured on this exact grid, that
+    construction gives +0.165 at (c+2s, c+2s) where the isotropic form gives -0.055.
+    """
+    som = make_som()
+    sigma = 3.0
+    h = som._mexicanhat((10, 10), sigma)
+    assert h[16, 16] < 0
+    assert h[16, 16] == pytest.approx(-0.054946916666, rel=1e-9)
+    assert h[19, 19] < 0
+
+
+def test_mexican_hat_has_no_positive_lobe_beyond_the_zero_crossing():
+    """Excitation is confined to the inner disc; everything beyond it inhibits or vanishes."""
+    som = make_som(x=41, y=41)
+    sigma = 4.0
+    h = som._mexicanhat((20, 20), sigma)
+    r = grid_radius((41, 41), (20, 20))
+    assert h[r > np.sqrt(2) * sigma + 1e-9].max() <= 0.0
+
+
+def test_mexican_hat_differs_from_the_separable_product():
+    """Guards against a regression to np.outer of two 1-D wavelets."""
+    som = make_som()
+    sigma = 3.0
+    dx = np.arange(21) - 10
+    ax = (1 - dx**2 / sigma**2) * np.exp(-(dx**2) / (2 * sigma**2))
+    separable = np.outer(ax, ax)
+    h = som._mexicanhat((10, 10), sigma)
+    assert not np.allclose(h, separable)
+    # The decisive difference: on the diagonal the separable product is excitatory where the
+    # isotropic form inhibits, because there both of its 1-D factors are negative.
+    assert separable[16, 16] > 0 > h[16, 16]
+    assert separable[16, 16] == pytest.approx(0.164840749998, rel=1e-9)
+    # The separable product is also anisotropic: (10, 14) and (13, 13) are not equidistant
+    # in its value despite both lying at grid radius 4 and 3*sqrt(2) respectively.
+    assert separable[10, 6] != pytest.approx(separable[13, 13])
+
+
+# --------------------------------------------------------------------------------------------
+# Bubble
+# --------------------------------------------------------------------------------------------
+
+
+def test_bubble_is_the_indicator_of_a_neighbourhood_set():
+    """Vrieze (1995) p. 85: N_i = {i' : d(i, i') <= rho}, the truncated inner lobe."""
+    som = make_som()
+    h = som._bubble((10, 10), 2.0)
+    assert set(np.unique(h)) <= {0, 1}
+    expected = np.zeros((21, 21), dtype=int)
+    expected[8:13, 8:13] = 1
+    np.testing.assert_array_equal(h, expected)
+
+
+# --------------------------------------------------------------------------------------------
+# Cyclic (toroidal) grids
+# --------------------------------------------------------------------------------------------
+
+
+@pytest.mark.parametrize("name", ["gaussian", "mexicanhat"])
+def test_cyclic_wraps_from_both_edges(name):
+    """Regression for the one-sided fold-back.
+
+    The original fold only handled dx > shape/2, leaving offsets below -shape/2 unfolded. A winner
+    at row 0 wrapped correctly (h[9] = 0.6065) while a winner at row 9 did not (h[0] = 0.0).
+    """
+    som = make_som(x=10, y=10, cyclic_x=True, cyclic_y=True, neighborhood_function=name)
+    low = som._neighborhood_function((0, 5), 1.0)[:, 5]
+    high = som._neighborhood_function((9, 5), 1.0)[:, 5]
+    assert low[9] == pytest.approx(low[1])
+    assert high[0] == pytest.approx(high[8])
+    np.testing.assert_allclose(high, np.roll(low, 9))
+
+
+@pytest.mark.parametrize("name", ["gaussian", "mexicanhat"])
+def test_cyclic_response_is_translation_invariant(name):
+    """On a torus every node is equivalent, so h depends only on the offset from the winner."""
+    som = make_som(x=12, y=12, cyclic_x=True, cyclic_y=True, neighborhood_function=name)
+    reference = som._neighborhood_function((0, 0), 2.0)
+    for cx in range(12):
+        for cy in range(12):
+            shifted = som._neighborhood_function((cx, cy), 2.0)
+            np.testing.assert_allclose(
+                shifted, np.roll(reference, (cx, cy), axis=(0, 1))
+            )
+
+
+def test_non_cyclic_grid_does_not_wrap():
+    som = make_som(x=10, y=10)
+    h = som._gaussian((0, 5), 1.0)[:, 5]
+    assert h[9] < 1e-12
+
+
+# --------------------------------------------------------------------------------------------
+# Selection and guards at the SOM level
+# --------------------------------------------------------------------------------------------
+
+
+def test_unknown_neighborhood_function_raises_value_error():
+    with pytest.raises(ValueError, match="neighborhood_function"):
+        python_som.SOM(x=8, y=8, input_len=3, neighborhood_function="mexican_hat")
+
+
+def test_batch_training_rejects_the_mexican_hat():
+    """Kohonen Eq. (8) divides by sum_j n_j h_ji, which is not sign-definite for a signed h."""
+    rng = np.random.default_rng(0)
+    data = rng.normal(size=(30, 3))
+    som = make_som(x=12, y=12, neighborhood_function="mexicanhat")
+    with pytest.raises(ValueError, match="batch"):
+        som.train(data, n_iteration=1, mode="batch")
+
+
+@pytest.mark.parametrize("mode", ["random", "sequential"])
+def test_stepwise_training_accepts_the_mexican_hat(mode):
+    rng = np.random.default_rng(0)
+    som = make_som(x=12, y=12, neighborhood_function="mexicanhat")
+    som.train(rng.normal(size=(30, 3)), n_iteration=30, mode=mode)
+    assert np.isfinite(som.get_weights()).all()