Merge pull request #80 from OVVO-Financial/claude/book-python-code-integration-xijzxt

Support unequal-length samples in stochastic dominance and nns_norm

Author: OVVO-Financial

docs/api_reference.md

@@ -315,7 +315,11 @@ NNS ANOVA-style comparison helper covering binary, multi-group, pairwise, and de
 
 Closest R API: `NNS.norm`.
 
-NNS normalization helper for numeric matrix-style inputs.
+NNS normalization helper. Accepts a 2-D matrix or, like R's list input, a
+sequence of 1-D vectors (one per variable). Equal-length vectors are
+column-stacked and normalized through the matrix path; unequal-length
+vectors force `linear=True` exactly as R does and return a list of scaled
+arrays.
 
 #### `nns_distance`
 
@@ -361,13 +365,17 @@ Rescales inputs using NNS conventions.
 
 Closest R APIs: `NNS.FSD`, `NNS.SSD`, and `NNS.TSD`.
 
-Compute first-, second-, and third-order stochastic dominance.
+Compute first-, second-, and third-order stochastic dominance. The two samples
+may differ in length; curves are evaluated on the merged threshold grid exactly
+as the R functions do.
 
 #### `fsd_uni`, `ssd_uni`, `tsd_uni`
 
 Closest R APIs: `NNS.FSD.uni`, `NNS.SSD.uni`, and `NNS.TSD.uni`.
 
-Univariate wrappers for stochastic dominance workflows.
+Univariate wrappers for stochastic dominance workflows. Unlike R's C++ `.uni`
+routines, unequal-length samples are supported with the same merged-grid
+semantics as the pairwise tests.
 
 #### `nns_sd_cluster`
 

docs/api_status.md

@@ -57,7 +57,7 @@ invariant, and property coverage.
 | Bootstrap/Monte Carlo: `nns_meboot`, `nns_mc` | implemented | medium | Deterministic diagnostics are parity-tested; exact stochastic replicate parity with R is not expected. |
 | Stochastic dominance/superiority: `fsd`, `ssd`, `tsd`, `.uni` wrappers, `nns_ss`, `nns_sd_cluster`, `sd_efficient_set` | implemented | medium | Public structures and deterministic paths are covered. SD uses exact pure-NumPy prefix-pair kernels plus a degree-1 discrete order-statistic matrix path; R's C++ core remains faster on full finance fixtures. Stochastic intervals use NNS Python RNG. |
 | ANOVA: `nns_anova` | implemented | high | Binary, multi-group, pairwise, and degenerate `NaN` conventions are covered. |
-| Normalization: `nns_norm` | implemented | high | Numeric matrix path is implemented. |
+| Normalization: `nns_norm` | implemented | high | Numeric matrix path and R's list-of-vectors path are implemented; unequal-length vectors force linear scaling as in R. |
 | Categorical helpers: `encode_factor_codes`, `factor_2_dummy`, `factor_2_dummy_fr`, `prepare_factor_predictors` | implemented | high | Explicit `levels=` / `factor_levels=` should be used to reproduce R factor ordering. `prepare_factor_predictors(...)` exposes the regression-ready full-rank design matrix path. |
 | Scalar differentiation: `nns_diff`, `dy_dx` | implemented | high | `dy_dx(..., eval_point="overall")` and numeric evaluation points are covered. |
 | Multivariate differentiation: `dy_d` | partial | medium-high | Scalar and vectorized point/distribution modes are covered on focused fixtures. Mixed derivatives are supported for two-regressor inputs where defined; multi-row matrix mixed derivatives use pointwise Python semantics rather than R's order-dependent list-matrix packing quirk. |

docs/conventions.md

@@ -73,7 +73,14 @@ NNS Python does not plot the dendrogram; it only returns the object data.
 The stochastic-dominance implementation is deliberately pure NumPy. It mirrors
 R's C++ SD core mathematically by sorting each column once, storing prefix sums,
 and evaluating dominance on each pair's merged threshold grid rather than on one
-global all-column grid. The full prefix-pair dominance matrix remains available
+global all-column grid. Pairwise tests (`fsd`, `ssd`, `tsd`, and the `*_uni`
+wrappers) accept samples of unequal length: each sample's curve is evaluated on
+the merged grid exactly as R's `NNS.FSD`/`NNS.SSD`/`NNS.TSD` compute
+`LPM(degree, sort(c(x, y)), sample)` per sample. R's C++ `.uni` walkers assume
+equal-length inputs, so for unequal lengths the Python `*_uni` wrappers are an
+intentional extension carrying the same merged-grid semantics. The
+matrix-based efficient-set and cluster routines operate on data columns and
+therefore remain equal-length by construction. The full prefix-pair dominance matrix remains available
 internally for verification and fallback. Large degree-1 discrete calls use an
 exact order-statistic dominance matrix: with equal-length empirical samples,
 one sample first-order stochastically dominates another exactly when every

src/nns/norm.py

@@ -1,16 +1,52 @@
 from __future__ import annotations
 
-from typing import cast
+from collections.abc import Sequence
+from typing import cast, overload
 
 import numpy as np
 from numpy.typing import NDArray
 
 from nns.dependence import nns_dep
 
 
-def nns_norm(x: NDArray[np.float64], linear: bool = False) -> NDArray[np.float64]:
-    """Normalize a numeric matrix following R's NNS.norm scaling."""
-    values = _as_matrix(x)
+@overload
+def nns_norm(x: NDArray[np.float64], linear: bool = ...) -> NDArray[np.float64]: ...
+
+
+@overload
+def nns_norm(
+    x: Sequence[NDArray[np.float64]],
+    linear: bool = ...,
+) -> NDArray[np.float64] | list[NDArray[np.float64]]: ...
+
+
+def nns_norm(
+    x: NDArray[np.float64] | Sequence[NDArray[np.float64]],
+    linear: bool = False,
+) -> NDArray[np.float64] | list[NDArray[np.float64]]:
+    """Normalize variables following R's NNS.norm scaling.
+
+    Two input conventions are supported, matching R's ``NNS.norm(X, ...)``:
+
+    * A 2-D array whose columns are variables. Returns the scaled 2-D array.
+    * A list or tuple of 1-D arrays, one per variable (R's list input; the
+      elements are variables/columns, not observation rows). Equal-length
+      vectors are column-stacked and normalized through the matrix path,
+      returning a 2-D array — mirroring R, where ``mapply`` simplifies the
+      equal-length list result to a matrix. Unequal-length vectors force
+      ``linear=True`` exactly as R does (dependence-based scale factors need
+      aligned columns) and return a list of scaled arrays, one per input
+      vector.
+    """
+    if isinstance(x, np.ndarray):
+        values = _as_matrix(x)
+    else:
+        series = [_as_vector(item, index) for index, item in enumerate(x)]
+        if not series:
+            raise ValueError("x must be non-empty.")
+        if len({item.size for item in series}) > 1:
+            return _norm_unequal_series(series)
+        values = _as_matrix(np.column_stack(series))
     means = np.mean(values, axis=0)
     means = means.copy()
     means[means == 0.0] = 1e-10
@@ -30,6 +66,14 @@ def nns_norm(x: NDArray[np.float64], linear: bool = False) -> NDArray[np.float64
     return scaled
 
 
+def _norm_unequal_series(series: list[NDArray[np.float64]]) -> list[NDArray[np.float64]]:
+    means = np.array([float(np.mean(item)) for item in series])
+    means[means == 0.0] = 1e-10
+    ratio_grid = means[:, np.newaxis] * (1.0 / means[np.newaxis, :])
+    scales = np.mean(ratio_grid, axis=0)
+    return [item * scale for item, scale in zip(series, scales, strict=True)]
+
+
 def _scale_factor(values: NDArray[np.float64]) -> NDArray[np.float64]:
     if values.shape[1] < 10:
         return cast(NDArray[np.float64], np.abs(np.corrcoef(values, rowvar=False)))
@@ -44,6 +88,17 @@ def _scale_factor(values: NDArray[np.float64]) -> NDArray[np.float64]:
     return deps
 
 
+def _as_vector(x: NDArray[np.float64], index: int) -> NDArray[np.float64]:
+    values = np.asarray(x, dtype=np.float64)
+    if values.ndim != 1:
+        raise ValueError(f"x[{index}] must be a 1D numeric vector.")
+    if values.size == 0:
+        raise ValueError(f"x[{index}] must be non-empty.")
+    if not np.all(np.isfinite(values)):
+        raise ValueError(f"x[{index}] must contain only finite values.")
+    return values
+
+
 def _as_matrix(x: NDArray[np.float64]) -> NDArray[np.float64]:
     values = np.asarray(x, dtype=np.float64)
     if values.ndim != 2:

src/nns/stochastic_dominance.py

@@ -4,6 +4,12 @@
 curve equality guard: equal LPM/CDF curves are non-dominance even when samples
 differ below meaningful double precision. Efficient-set output follows the R
 C++ routine's LPM-at-global-maximum ordering and original-index tie break.
+
+Pairwise tests accept samples of unequal length: each sample's LPM/CDF curve is
+evaluated on the merged threshold grid, exactly as R's NNS.FSD/NNS.SSD/NNS.TSD
+compute ``LPM(degree, sort(c(x, y)), sample)`` per sample. (R's C++ ``.uni``
+walkers assume equal lengths; the Python ``*_uni`` wrappers extend the same
+merged-grid semantics to unequal lengths.)
 """
 
 from __future__ import annotations
@@ -55,7 +61,7 @@ class _SDOrderStatPrecomputed:
 
 
 def fsd(x: NDArray[np.float64], y: NDArray[np.float64]) -> int:
-    """First-order stochastic dominance."""
+    """First-order stochastic dominance; ``x`` and ``y`` may differ in length."""
     x_values = _as_sd_values(x, "x")
     y_values = _as_sd_values(y, "y")
     return _sd_result(x_values, y_values, 1)
@@ -70,7 +76,7 @@ def fsd_uni(x: NDArray[np.float64], y: NDArray[np.float64], type: str = "discret
 
 
 def ssd(x: NDArray[np.float64], y: NDArray[np.float64]) -> int:
-    """Second-order stochastic dominance."""
+    """Second-order stochastic dominance; ``x`` and ``y`` may differ in length."""
     x_values = _as_sd_values(x, "x")
     y_values = _as_sd_values(y, "y")
     return _sd_result(x_values, y_values, 2)
@@ -84,7 +90,7 @@ def ssd_uni(x: NDArray[np.float64], y: NDArray[np.float64]) -> int:
 
 
 def tsd(x: NDArray[np.float64], y: NDArray[np.float64]) -> int:
-    """Third-order stochastic dominance."""
+    """Third-order stochastic dominance; ``x`` and ``y`` may differ in length."""
     x_values = _as_sd_values(x, "x")
     y_values = _as_sd_values(y, "y")
     return _sd_result(x_values, y_values, 3)
@@ -754,8 +760,6 @@ def _dominates_uni(
     *,
     discrete: bool,
 ) -> bool:
-    if x.size != y.size:
-        raise ValueError("x and y must have the same length.")
     if np.array_equal(np.sort(x), np.sort(y)):
         return False
     if np.min(x) < np.min(y):

tests/invariants/test_norm.py

@@ -34,3 +34,60 @@ def test_nonlinear_nns_norm_preserves_shape_for_wide_matrix() -> None:
 
     assert result.shape == x.shape
     assert np.all(np.isfinite(result))
+
+
+def test_equal_length_sequence_matches_matrix_path() -> None:
+    x = np.linspace(1.0, 3.0, 50)
+    y = np.linspace(2.0, 8.0, 50)
+    z = np.linspace(10.0, 20.0, 50)
+    matrix = np.column_stack((x, y, z))
+
+    for linear in (False, True):
+        from_sequence = nns_norm([x, y, z], linear=linear)
+        from_matrix = nns_norm(matrix, linear=linear)
+        assert isinstance(from_sequence, np.ndarray)
+        np.testing.assert_allclose(from_sequence, from_matrix)
+
+
+def test_unequal_length_sequence_forces_linear_scaling() -> None:
+    vec1 = np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0])
+    vec2 = np.array([10.0, 20.0, 30.0, 40.0, 50.0, 60.0])
+    vec3 = np.array([0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3])
+
+    result = nns_norm([vec1, vec2, vec3])
+
+    assert isinstance(result, list)
+    assert [item.size for item in result] == [7, 6, 9]
+
+    # Linear scaling equalizes every scaled mean at the grand mean of means,
+    # and the linear flag is forced regardless of its passed value (as in R).
+    grand_mean = np.mean([vec1.mean(), vec2.mean(), vec3.mean()])
+    for item in result:
+        np.testing.assert_allclose(item.mean(), grand_mean)
+    forced = nns_norm([vec1, vec2, vec3], linear=True)
+    for got, expected in zip(result, forced, strict=True):
+        np.testing.assert_allclose(got, expected)
+
+
+def test_zero_sum_length_differences_detected_as_unequal() -> None:
+    # Lengths (5, 6, 5) have pairwise diffs summing to zero; they must still
+    # take the unequal-length path rather than being treated as a matrix.
+    vec1 = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
+    vec2 = np.array([10.0, 20.0, 30.0, 40.0, 50.0, 60.0])
+    vec3 = np.array([0.5, 0.6, 0.7, 0.8, 0.9])
+
+    result = nns_norm([vec1, vec2, vec3])
+
+    assert isinstance(result, list)
+    assert [item.size for item in result] == [5, 6, 5]
+
+
+def test_sequence_input_validation() -> None:
+    import pytest
+
+    with pytest.raises(ValueError, match="non-empty"):
+        nns_norm([])
+    with pytest.raises(ValueError, match=r"x\[1\] must be a 1D"):
+        nns_norm([np.ones(3), np.ones((3, 2))])
+    with pytest.raises(ValueError, match=r"x\[0\] must contain only finite"):
+        nns_norm([np.array([1.0, np.nan]), np.ones(3)])

tests/invariants/test_stochastic_dominance.py

@@ -1,8 +1,32 @@
 from __future__ import annotations
 
 import numpy as np
+import pytest
+from numpy.typing import NDArray
 
-from nns import fsd, ssd, tsd
+from nns import fsd, fsd_uni, lpm, ssd, ssd_uni, tsd, tsd_uni
+
+
+def _r_sd_reference(x: NDArray[np.float64], y: NDArray[np.float64], degree: int) -> int:
+    """R's NNS.FSD/NNS.SSD/NNS.TSD decision rule, transcribed literally."""
+    grid = np.sort(np.concatenate((x, y)))
+    if degree == 1:
+        # LPM.ratio(0, grid, sample) == LPM(0, grid, sample) == ECDF
+        curve_x = np.asarray(lpm(0, grid, x), dtype=np.float64)
+        curve_y = np.asarray(lpm(0, grid, y), dtype=np.float64)
+    else:
+        curve_x = np.asarray(lpm(degree - 1, grid, x), dtype=np.float64)
+        curve_y = np.asarray(lpm(degree - 1, grid, y), dtype=np.float64)
+
+    mean_ok_xy = degree == 1 or float(np.mean(x)) >= float(np.mean(y))
+    mean_ok_yx = degree == 1 or float(np.mean(y)) >= float(np.mean(x))
+    curves_identical = np.array_equal(curve_x, curve_y)
+
+    if not np.any(curve_x > curve_y) and x.min() >= y.min() and mean_ok_xy and not curves_identical:
+        return 1
+    if not np.any(curve_y > curve_x) and y.min() >= x.min() and mean_ok_yx and not curves_identical:
+        return -1
+    return 0
 
 
 def test_sd_antisymmetry() -> None:
@@ -29,3 +53,60 @@ def test_self_does_not_dominate() -> None:
     assert fsd(x, x) == 0
     assert ssd(x, x) == 0
     assert tsd(x, x) == 0
+
+
+def test_unequal_length_shifted_samples_dominate() -> None:
+    x = np.array([2.0, 3.0, 4.0, 5.0, 6.0])
+    y = np.array([1.0, 2.0, 3.0])
+
+    assert fsd(x, y) == 1
+    assert ssd(x, y) == 1
+    assert tsd(x, y) == 1
+    assert fsd(y, x) == -1
+
+    assert fsd_uni(x, y) == 1
+    assert fsd_uni(x, y, "continuous") == 1
+    assert ssd_uni(x, y) == 1
+    assert tsd_uni(x, y) == 1
+    assert fsd_uni(y, x) == 0
+    assert ssd_uni(y, x) == 0
+    assert tsd_uni(y, x) == 0
+
+
+def test_unequal_length_antisymmetry() -> None:
+    rng = np.random.default_rng(7)
+    x = rng.normal(1.0, 1.0, 37)
+    y = rng.normal(0.0, 1.0, 61)
+
+    assert fsd(x, y) == -fsd(y, x)
+    assert ssd(x, y) == -ssd(y, x)
+    assert tsd(x, y) == -tsd(y, x)
+
+
+def test_mtcars_transmission_groups_match_r() -> None:
+    # mtcars mpg split by transmission: R's NNS.FSD returns "X FSD Y" for
+    # (manual, auto) despite the samples having different lengths (13 vs 19).
+    auto_mpg = np.array(
+        [21.4, 18.7, 18.1, 14.3, 24.4, 22.8, 19.2, 17.8, 16.4, 17.3,
+         15.2, 10.4, 10.4, 14.7, 21.5, 15.5, 15.2, 13.3, 19.2]
+    )
+    manual_mpg = np.array(
+        [21.0, 21.0, 22.8, 32.4, 30.4, 33.9, 27.3, 26.0, 30.4, 15.8, 19.7, 15.0, 21.4]
+    )
+
+    assert fsd(manual_mpg, auto_mpg) == 1
+    assert fsd(auto_mpg, manual_mpg) == -1
+    assert fsd_uni(manual_mpg, auto_mpg) == 1
+
+
+@pytest.mark.parametrize("degree", [1, 2, 3])
+@pytest.mark.parametrize("seed", range(8))
+def test_unequal_length_matches_r_decision_rule(degree: int, seed: int) -> None:
+    rng = np.random.default_rng(seed)
+    sizes = rng.integers(3, 60, size=2)
+    shift = rng.uniform(-0.5, 0.5)
+    x = rng.normal(shift, 1.0, int(sizes[0]))
+    y = rng.normal(0.0, rng.uniform(0.5, 1.5), int(sizes[1]))
+
+    function = {1: fsd, 2: ssd, 3: tsd}[degree]
+    assert function(x, y) == _r_sd_reference(x, y, degree)