ci(d0): run Tests/eval non-vacuously, and fix what that surfaced

.github/workflows/ci.yml

@@ -79,3 +79,82 @@ jobs:
         # collision/flake risk on a headless runner. The OSC *decode* logic is
         # covered socket-free by MuseOSCDecoderTests / MindMonitorDecoderTests,
         # which still run. Run the socket integration test locally.
+
+  # Tests/eval had no CI at all before this job. Seven modules, and only
+  # test_adr_references.py ran (via the direct python3 step above) — the other
+  # six executed nowhere, which is how three failing assertions and an
+  # import-time crash survived in a green repository.
+  #
+  # pytest, not unittest, and the choice is forced rather than stylistic:
+  # six of the seven modules are bare `def test_*` with no `unittest.TestCase`,
+  # so `python -m unittest` collects ZERO from them, and `unittest discover`
+  # fails outright here (`Tests/eval` has no __init__.py). Either way the suite
+  # would report success while running nothing.
+  python-eval:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+
+      - uses: actions/setup-python@v5
+        with:
+          python-version: "3.11"
+
+      - name: Install Python dependencies
+        run: |
+          # scipy is required, not optional: three modules import it
+          # transitively and fail collection without it. pandas is a harness
+          # dependency of the session/channel-quality frames.
+          #
+          # Majors are bounded because this gate's semantics depend on them:
+          # the trapz/trapezoid fix is numpy-2 specific, and the Mann-Whitney
+          # contract depends on SciPy's exact/asymptotic behaviour. An
+          # unconstrained install would let a major bump change a p-value
+          # without any commit here.
+          pip install \
+            'pytest>=8,<9' \
+            'numpy>=2,<3' \
+            'scipy>=1.14,<2' \
+            'pandas>=2,<3'
+
+      # A --filter or --ignore that matches nothing exits 0, and a collection
+      # error can silently shrink the suite. Assert the floor before running so
+      # a module that stops being collected fails loudly instead of vanishing.
+      - name: Assert Tests/eval is discoverable
+        run: |
+          python -m pytest Tests/eval --collect-only -q \
+            --ignore=Tests/eval/test_joint_embedding.py > /tmp/collect.txt 2>&1 || true
+          tail -1 /tmp/collect.txt
+          collected=$(grep -oE '^[0-9]+ tests? collected' /tmp/collect.txt | grep -oE '^[0-9]+' | head -1)
+          errors=$(grep -cE '^ERROR ' /tmp/collect.txt || true)
+          echo "discovered: Tests/eval=${collected:-0} collection_errors=${errors:-0}"
+          [ "${errors:-0}" -eq 0 ] \
+            || { echo "::error::${errors} module(s) failed to collect"; exit 1; }
+          # Floor sits below the current 63 so ordinary additions never trip it.
+          [ "${collected:-0}" -ge 55 ] \
+            || { echo "::error::Tests/eval collected ${collected:-0}, expected >= 55"; exit 1; }
+
+      - name: Test (eval suites)
+        run: |
+          # test_joint_embedding.py imports mlx (Apple Silicon only) and needs a
+          # converted model; it is the one module excluded, by name and reason.
+          python -m pytest Tests/eval -v \
+            --ignore=Tests/eval/test_joint_embedding.py
+
+      # NeuralComposeEEG is not on main yet — it arrives with lane B. Absent is
+      # expected and must be said out loud; present-but-empty is a defect, since
+      # a package whose tests silently collect nothing is the failure this whole
+      # job exists to prevent.
+      - name: Assert NeuralComposeEEG is absent or non-vacuous
+        run: |
+          if [ ! -d NeuralComposeEEG ]; then
+            echo "NeuralComposeEEG/ absent on this branch — expected until lane B integrates it"
+            exit 0
+          fi
+          echo "NeuralComposeEEG/ present — its suite must not be vacuous"
+          n=$(PYTHONPATH=NeuralComposeEEG/src python -m pytest NeuralComposeEEG/tests \
+                --collect-only -q 2>/dev/null \
+                | grep -oE '^[0-9]+ tests? collected' | grep -oE '^[0-9]+' | head -1)
+          echo "discovered: NeuralComposeEEG=${n:-0}"
+          [ "${n:-0}" -gt 0 ] \
+            || { echo "::error::NeuralComposeEEG/ is present but collected 0 tests"; exit 1; }
+          PYTHONPATH=NeuralComposeEEG/src python -m pytest NeuralComposeEEG/tests -v

Evaluation/scripts/analyze_results.py

@@ -44,10 +44,31 @@ def cohens_d(a, b):
     a, b = np.array(a, dtype=float), np.array(b, dtype=float)
     if len(a) < 2 or len(b) < 2:
         return float("nan")
-    pooled_std = math.sqrt((a.var(ddof=1) + b.var(ddof=1)) / 2)
+    # Degrees-of-freedom-weighted pooled variance. The unweighted
+    # (var_a + var_b) / 2 form coincides with this only when the groups are the
+    # same size, and this function accepts arbitrary lengths — so every
+    # unequal-n caller was dividing by the wrong denominator.
+    n_a, n_b = len(a), len(b)
+    pooled_variance = (
+        (n_a - 1) * a.var(ddof=1) + (n_b - 1) * b.var(ddof=1)
+    ) / (n_a + n_b - 2)
+    pooled_std = math.sqrt(pooled_variance)
+    mean_difference = a.mean() - b.mean()
     if pooled_std == 0:
-        return 0.0
-    return float((a.mean() - b.mean()) / pooled_std)
+        # PROJECT POLICY, not a universal definition. With both samples constant
+        # the denominator is zero and Cohen's d is strictly undefined. We report
+        # the limiting value — signed infinity — because 0.0 claimed "no effect"
+        # for the most extreme separation representable, and nan would discard
+        # the direction. Identical constants are the one genuinely zero case.
+        #
+        # Callers that PERSIST this must convert at the boundary: JSON has no
+        # Infinity literal (`json.dumps` emits a bare `Infinity`, which strict
+        # parsers reject) and CSV has no convention at all. Decide on a
+        # representation before writing it, rather than assuming the reader copes.
+        if mean_difference == 0:
+            return 0.0
+        return float("inf") if mean_difference > 0 else float("-inf")
+    return float(mean_difference / pooled_std)
 
 
 def safe_mean(values):

Evaluation/scripts/embedding_analyze.py

@@ -54,19 +54,56 @@ def cohens_d(a, b):
     a, b = np.array(a, dtype=float), np.array(b, dtype=float)
     if len(a) < 2 or len(b) < 2:
         return float("nan")
-    pooled_std = math.sqrt((a.var(ddof=1) + b.var(ddof=1)) / 2)
+    # Degrees-of-freedom-weighted pooled variance. The unweighted
+    # (var_a + var_b) / 2 form coincides with this only when the groups are the
+    # same size, and this function accepts arbitrary lengths — so every
+    # unequal-n caller was dividing by the wrong denominator.
+    n_a, n_b = len(a), len(b)
+    pooled_variance = (
+        (n_a - 1) * a.var(ddof=1) + (n_b - 1) * b.var(ddof=1)
+    ) / (n_a + n_b - 2)
+    pooled_std = math.sqrt(pooled_variance)
+    mean_difference = a.mean() - b.mean()
     if pooled_std == 0:
-        return 0.0
-    return float((a.mean() - b.mean()) / pooled_std)
-
+        # PROJECT POLICY, not a universal definition. With both samples constant
+        # the denominator is zero and Cohen's d is strictly undefined. We report
+        # the limiting value — signed infinity — because 0.0 claimed "no effect"
+        # for the most extreme separation representable, and nan would discard
+        # the direction. Identical constants are the one genuinely zero case.
+        #
+        # Callers that PERSIST this must convert at the boundary: JSON has no
+        # Infinity literal (`json.dumps` emits a bare `Infinity`, which strict
+        # parsers reject) and CSV has no convention at all. Decide on a
+        # representation before writing it, rather than assuming the reader copes.
+        if mean_difference == 0:
+            return 0.0
+        return float("inf") if mean_difference > 0 else float("-inf")
+    return float(mean_difference / pooled_std)
+
+
+# Above this per-group size the exact Mann-Whitney null distribution is too
+# expensive to enumerate, so the asymptotic approximation is used. Stated here
+# rather than inherited from SciPy's `auto` threshold, which is a library
+# implementation detail and not a contract.
+_EXACT_MAX_N = 20
 
 def mann_whitney_u(a, b):
     a, b = np.array(a, dtype=float), np.array(b, dtype=float)
     if len(a) < 2 or len(b) < 2:
         return {"u_statistic": float("nan"), "p_value": float("nan"),
                 "effect_size_r": float("nan"), "n_a": len(a), "n_b": len(b)}
     try:
-        u_stat, p_value = sp_stats.mannwhitneyu(a, b, alternative="two-sided")
+        # Method chosen explicitly rather than inherited from SciPy's `auto`,
+        # whose selection rule has changed across versions and would silently
+        # move a reported p-value between releases. Exact is only valid without
+        # ties; it is also combinatorial, so it is capped by sample size and the
+        # cap is a stated policy rather than a SciPy default.
+        combined = np.concatenate([a, b])
+        has_ties = len(np.unique(combined)) < len(combined)
+        too_large = max(len(a), len(b)) > _EXACT_MAX_N
+        method = "asymptotic" if (has_ties or too_large) else "exact"
+        u_stat, p_value = sp_stats.mannwhitneyu(
+            a, b, alternative="two-sided", method=method)
         n = len(a) + len(b)
         r = 1 - (2 * u_stat) / (n * (n - 1) / 2) if n > 1 else float("nan")
         return {

Evaluation/scripts/eval_stats.py

@@ -36,11 +36,38 @@ def cohens_d(a, b):
     a, b = np.array(a, dtype=float), np.array(b, dtype=float)
     if len(a) < 2 or len(b) < 2:
         return float("nan")
-    pooled_std = math.sqrt((a.var(ddof=1) + b.var(ddof=1)) / 2)
+    # Degrees-of-freedom-weighted pooled variance. The unweighted
+    # (var_a + var_b) / 2 form coincides with this only when the groups are the
+    # same size, and this function accepts arbitrary lengths — so every
+    # unequal-n caller was dividing by the wrong denominator.
+    n_a, n_b = len(a), len(b)
+    pooled_variance = (
+        (n_a - 1) * a.var(ddof=1) + (n_b - 1) * b.var(ddof=1)
+    ) / (n_a + n_b - 2)
+    pooled_std = math.sqrt(pooled_variance)
+    mean_difference = a.mean() - b.mean()
     if pooled_std == 0:
-        return 0.0
-    return float((a.mean() - b.mean()) / pooled_std)
-
+        # PROJECT POLICY, not a universal definition. With both samples constant
+        # the denominator is zero and Cohen's d is strictly undefined. We report
+        # the limiting value — signed infinity — because 0.0 claimed "no effect"
+        # for the most extreme separation representable, and nan would discard
+        # the direction. Identical constants are the one genuinely zero case.
+        #
+        # Callers that PERSIST this must convert at the boundary: JSON has no
+        # Infinity literal (`json.dumps` emits a bare `Infinity`, which strict
+        # parsers reject) and CSV has no convention at all. Decide on a
+        # representation before writing it, rather than assuming the reader copes.
+        if mean_difference == 0:
+            return 0.0
+        return float("inf") if mean_difference > 0 else float("-inf")
+    return float(mean_difference / pooled_std)
+
+
+# Above this per-group size the exact Mann-Whitney null distribution is too
+# expensive to enumerate, so the asymptotic approximation is used. Stated here
+# rather than inherited from SciPy's `auto` threshold, which is a library
+# implementation detail and not a contract.
+_EXACT_MAX_N = 20
 
 def mann_whitney_u(a, b):
     """Mann-Whitney U test."""
@@ -49,7 +76,17 @@ def mann_whitney_u(a, b):
         return {"u_statistic": float("nan"), "p_value": float("nan"),
                 "effect_size_r": float("nan"), "n_a": len(a), "n_b": len(b)}
     try:
-        u_stat, p_value = sp_stats.mannwhitneyu(a, b, alternative="two-sided")
+        # Method chosen explicitly rather than inherited from SciPy's `auto`,
+        # whose selection rule has changed across versions and would silently
+        # move a reported p-value between releases. Exact is only valid without
+        # ties; it is also combinatorial, so it is capped by sample size and the
+        # cap is a stated policy rather than a SciPy default.
+        combined = np.concatenate([a, b])
+        has_ties = len(np.unique(combined)) < len(combined)
+        too_large = max(len(a), len(b)) > _EXACT_MAX_N
+        method = "asymptotic" if (has_ties or too_large) else "exact"
+        u_stat, p_value = sp_stats.mannwhitneyu(
+            a, b, alternative="two-sided", method=method)
         n = len(a) + len(b)
         r = 1 - (2 * u_stat) / (n * (n - 1) / 2) if n > 1 else float("nan")
         return {

Evaluation/scripts/statistical_analysis.py

@@ -60,11 +60,38 @@ def cohens_d(a, b):
     a, b = np.array(a, dtype=float), np.array(b, dtype=float)
     if len(a) < 2 or len(b) < 2:
         return float("nan")
-    pooled_std = math.sqrt((a.var(ddof=1) + b.var(ddof=1)) / 2)
+    # Degrees-of-freedom-weighted pooled variance. The unweighted
+    # (var_a + var_b) / 2 form coincides with this only when the groups are the
+    # same size, and this function accepts arbitrary lengths — so every
+    # unequal-n caller was dividing by the wrong denominator.
+    n_a, n_b = len(a), len(b)
+    pooled_variance = (
+        (n_a - 1) * a.var(ddof=1) + (n_b - 1) * b.var(ddof=1)
+    ) / (n_a + n_b - 2)
+    pooled_std = math.sqrt(pooled_variance)
+    mean_difference = a.mean() - b.mean()
     if pooled_std == 0:
-        return 0.0
-    return float((a.mean() - b.mean()) / pooled_std)
-
+        # PROJECT POLICY, not a universal definition. With both samples constant
+        # the denominator is zero and Cohen's d is strictly undefined. We report
+        # the limiting value — signed infinity — because 0.0 claimed "no effect"
+        # for the most extreme separation representable, and nan would discard
+        # the direction. Identical constants are the one genuinely zero case.
+        #
+        # Callers that PERSIST this must convert at the boundary: JSON has no
+        # Infinity literal (`json.dumps` emits a bare `Infinity`, which strict
+        # parsers reject) and CSV has no convention at all. Decide on a
+        # representation before writing it, rather than assuming the reader copes.
+        if mean_difference == 0:
+            return 0.0
+        return float("inf") if mean_difference > 0 else float("-inf")
+    return float(mean_difference / pooled_std)
+
+
+# Above this per-group size the exact Mann-Whitney null distribution is too
+# expensive to enumerate, so the asymptotic approximation is used. Stated here
+# rather than inherited from SciPy's `auto` threshold, which is a library
+# implementation detail and not a contract.
+_EXACT_MAX_N = 20
 
 def mann_whitney_u(a, b):
     """Mann-Whitney U test — non-parametric, no normality assumption."""
@@ -73,7 +100,17 @@ def mann_whitney_u(a, b):
         return {"u_statistic": float("nan"), "p_value": float("nan"),
                 "effect_size_r": float("nan"), "n_a": len(a), "n_b": len(b)}
     try:
-        u_stat, p_value = sp_stats.mannwhitneyu(a, b, alternative="two-sided")
+        # Method chosen explicitly rather than inherited from SciPy's `auto`,
+        # whose selection rule has changed across versions and would silently
+        # move a reported p-value between releases. Exact is only valid without
+        # ties; it is also combinatorial, so it is capped by sample size and the
+        # cap is a stated policy rather than a SciPy default.
+        combined = np.concatenate([a, b])
+        has_ties = len(np.unique(combined)) < len(combined)
+        too_large = max(len(a), len(b)) > _EXACT_MAX_N
+        method = "asymptotic" if (has_ties or too_large) else "exact"
+        u_stat, p_value = sp_stats.mannwhitneyu(
+            a, b, alternative="two-sided", method=method)
         # Rank-biserial correlation effect size
         n = len(a) + len(b)
         r = 1 - (2 * u_stat) / (n * (n - 1) / 2) if n > 1 else float("nan")

Scripts/eeg_spectral.py

@@ -33,10 +33,29 @@
 }
 
 _EPS = 1e-8
-# np.trapezoid is the numpy>=2.0 name; np.trapz is the numpy<2.0 name. This
-# environment pins numpy 1.26 (has only trapz); support both so the module is
-# portable if the calibration venv is later bumped to numpy 2.x.
-_trapz = getattr(np, "trapezoid", getattr(np, "trapz"))
+
+
+def _resolve_trapezoid(namespace):
+    """Pick the integrator NumPy actually exposes.
+
+    `trapezoid` is the numpy>=2.0 name and `trapz` the numpy<2.0 one. Resolving
+    it lazily is the whole point: the previous form was
+
+        getattr(np, "trapezoid", getattr(np, "trapz"))
+
+    and Python evaluates that default *before* the lookup it is a fallback for.
+    On numpy 2.x, which removed `trapz` outright, it raised `AttributeError` at
+    import — on precisely the version the fallback existed to support. Every
+    `Tests/eval` module reaching this file failed collection, which is how it
+    survived: nothing in CI imported it.
+    """
+    trapezoid = getattr(namespace, "trapezoid", None)
+    if trapezoid is not None:
+        return trapezoid
+    return namespace.trapz
+
+
+_trapz = _resolve_trapezoid(np)
 
 
 def band_power(freqs: np.ndarray, psd: np.ndarray, band: tuple[float, float]) -> float:

Tests/eval/test_embedding_space.py

@@ -21,11 +21,35 @@ def test_cka_identical():
 
 
 def test_cka_independent():
+    """Independent data does NOT give CKA near zero when d >> n.
+
+    This is the biased linear CKA estimator (Kornblith et al.). With 50 samples
+    against 128 and 256 features it inflates severely — ~0.78 here — because the
+    Gram matrices of independent high-dimensional Gaussians are far from
+    orthogonal at this sample size. The previous `< 0.3` described an unbiased
+    estimator this code does not implement, and the test never ran to contradict
+    it.
+
+    Pinned as an upper bound rather than an equality: the point is that the
+    inflation is real and bounded, so a future switch to an unbiased estimator
+    (which would drop this toward 0) fails here and forces the expectation to be
+    revisited deliberately.
+    """
     np.random.seed(42)
     X = np.random.randn(50, 128)
     Y = np.random.randn(50, 256)
     score = cka(X, Y)
-    assert abs(score) < 0.3
+    assert 0.0 <= score <= 1.0
+    assert score > 0.5, "biased CKA is expected to inflate at d >> n"
+
+
+def test_cka_inflation_shrinks_as_samples_grow():
+    """The inflation is a sample-size artefact, not a property of the data:
+    raising n against fixed d moves the independent-data score down."""
+    np.random.seed(42)
+    few = cka(np.random.randn(50, 128), np.random.randn(50, 128))
+    many = cka(np.random.randn(400, 128), np.random.randn(400, 128))
+    assert many < few
 
 
 def test_svcca_identical():