Fresh R NNS 13.0 parity regeneration and corrective Python fixes

.gitignore

@@ -14,5 +14,12 @@ __pycache__/
 # Virtual environments
 .venv/
 
-# R cache lock
+# R cache lock and fresh-regeneration backup
 tests/_r_cache.lock
+tests/_r_cache.json.bak
+
+# Vendored R NNS local-install build artifacts
+tools/NNS/src/*.o
+
+# Stray R plotting artifact from headless regeneration
+Rplots.pdf

docs/r13_cache_regeneration.md

@@ -0,0 +1,170 @@
+# R NNS 13.0 parity cache regeneration provenance
+
+This document records the fresh, from-empty regeneration of
+`tests/_r_cache.json` against live vendored R NNS 13.0.
+
+## Environment
+
+| Item | Value |
+| --- | --- |
+| Date of fresh regeneration | 2026-06-13 (UTC) |
+| OS | Ubuntu 24.04.4 LTS (Linux 6.18.5 x86_64) |
+| Python | 3.11.15 |
+| R | R version 4.3.3 (2024-02-29) "Angel Food Cake" |
+| Rscript | Rscript (R) version 4.3.3 (2024-02-29) |
+| Vendored NNS source | `tools/NNS` (extracted package directory, preferred over `tools/NNS_13.0.tar.gz`) |
+| `packageVersion("NNS")` | `13.0` |
+
+R NNS itself was installed exclusively from the vendored repository source via
+`scripts/install_local_r_nns.py` (`R CMD INSTALL tools/NNS`), never from CRAN.
+R package dependencies required to load vendored NNS (`data.table`,
+`doParallel`, `foreach`, `Rcpp`, `RcppParallel`, `rgl`, `xts`, `zoo`,
+`jsonlite`) were installed as Ubuntu binary packages; `Rfast` (plus its `zigg`
+dependency) was built from the upstream GitHub release source
+`RfastOfficial/Rfast` tag `v2.1.5.1-apollo` because CRAN was unreachable in
+the regeneration environment. Only dependencies came from external archives;
+NNS came from `tools/NNS`.
+
+## Commands run
+
+```bash
+python scripts/install_local_r_nns.py
+Rscript -e "suppressPackageStartupMessages(library(NNS)); cat(as.character(packageVersion('NNS')))"
+# printed: 13.0
+
+# Fresh regeneration: moves tests/_r_cache.json to tests/_r_cache.json.bak,
+# starts from an EMPTY cache, and repopulates every entry with a live R call.
+python scripts/regenerate_r_cache.py --fresh -- -n 0 tests/parity
+
+# Replay checks after regeneration
+PYNNS_R_CACHE_ONLY=1 python -m pytest -q -n 0 tests/parity
+python -m pytest -q tests/parity/test_r13_smoke.py
+python -m pytest -q tests/invariants
+ruff check .
+mypy
+python -m build
+```
+
+All cache-only/offline toggles (`PYNNS_R_CACHE_ONLY`, `NNS_R_CACHE_ONLY`,
+`PYNNS_OFFLINE`, `NNS_OFFLINE`, `CI`) were unset for the regeneration run.
+The `--fresh` mode added to `scripts/regenerate_r_cache.py` in this change
+refuses to run in CI, verifies a live local R NNS 13.0 install before touching
+anything, then moves the existing cache aside so no existing entry can be
+reused during regeneration.
+
+The fresh regeneration was run twice: once to surface the true Python parity
+gaps against live R 13.0 (216 test failures), and a second time from an empty
+cache after the Python fixes below, confirming every parity test passes
+against entries produced exclusively by live R calls.
+
+## Results
+
+| Item | Value |
+| --- | --- |
+| `nns_version` | `13.0` |
+| `schema_version` | `1` |
+| Final entry count | 2385 |
+| Prior entry count | 2406 |
+
+The entry count dropped from 2406 to 2385 because 21 entries in the old
+committed cache are stale: an instrumented cache-only replay confirmed that no
+current parity test computes those 21 keys, so a from-empty regeneration never
+recreates them. All 2385 keys requested by the current test suite were
+regenerated from live R.
+
+### Deterministic cache entries that changed
+
+456 of the 2384 shared keys changed value relative to the previously committed
+cache, summarized by NNS function (mapped by instrumenting the cache-key
+computation during a full replay):
+
+| Function | Changed entries |
+| --- | --- |
+| `PM.matrix` | 216 |
+| `NNS.reg` | 147 |
+| `NNS.boost` (numeric harness) | 25 |
+| `LPM.VaR` | 12 |
+| `UPM.VaR` | 12 |
+| `NNS.ANOVA` (custom harness) | 12 |
+| `NNS.M.reg` | 10 |
+| `NNS.stack` (numeric harness) | 6 |
+| `NNS.dep` | 3 |
+| `NNS.caus` | 3 |
+| `NNS.moments` | 2 |
+| `dy.d` (scalar harnesses) | 4 |
+| factor-predictor harnesses (boost/stack/reg) | 4 |
+
+The fresh live R 13.0 values were treated as authoritative in every case.
+
+### Python parity fixes required
+
+Fresh regeneration surfaced 216 failing parity tests. All of them traced to
+four deterministic divergences between NNS Python and vendored R NNS 13.0, and
+Python was fixed to match R in each case:
+
+1. **`LPM.VaR` / `UPM.VaR` integer degrees 1–4** (`src/nns/var.py`):
+   R 13.0 replaced the `optimize()` search with an exact polynomial
+   root-finding inversion on the located order-statistic interval
+   (`.NNS_LPM_VaR_integer`). Python now ports that algorithm (prefix power
+   sums, break-ratio interval location, Brent root with
+   `tol = .Machine$double.eps^0.5`, and the same boundary/fallback handling).
+   This also fixed the `NNS.reg`/`NNS.M.reg` confidence intervals and
+   stack/boost prediction intervals built on these helpers.
+2. **Distance-kernel lognormal rank weight** (`src/nns/distance.py`):
+   R 13.0's bulk path kernels (`NNS_distance_path_single_parallel_cpp`, used
+   for fitted values when `n.best > 1` and for multi-point estimates) use the
+   population sd of ranks `sqrt((k^2 - 1) / 12)` in the lognormal weight,
+   while the single-point `NNS_distance_cpp` kernel keeps the sample sd.
+   Python now provides `nns_distance_path_single_bulk` mirroring the bulk
+   kernel and keeps `nns_distance` on the single-point formulation.
+3. **`NNS.M.reg` out-of-hull multi-point extrapolation**
+   (`src/nns/multivariate_regression.py`): R 13.0 vectorized the multi-point
+   outsider path (bulk kernel estimates, `pmax(distance, 1e-10)` gradient
+   guards) and removed the old dims-dropping behavior for a single outsider
+   row that Python previously emulated. Python now mirrors the new path.
+4. **`NNS.ARMA` numeric multi-lag seasonal weighting** (`src/nns/arma.py`):
+   R weights each numeric seasonal factor by reversing the series with the
+   factor's *position* in the `seasonal.factor` vector
+   (`variable[seq(length(variable), 1, -i)]`), not its lag value. Python now
+   matches, which resolved the two previously `xfail`-ed Sunspots ARMA and
+   macro-like VAR practical examples; both now pass against live R and the
+   `xfail` markers were removed.
+
+Structural test updates justified by fresh R output: R 13.0's `NNS.boost`
+returns only `results`, `pred.int`, `feature.weights`, and
+`feature.frequency` (no `n.best`), so the Python return dictionary and the
+test assertions that expected `n.best` were updated, and the NaN handling of
+final boost estimates now matches R (`stack` falls back to `reg` only when
+absent; NaNs are filled with the gravity of the remaining estimates).
+
+The only remaining `xfail` is the balanced Iris `NNS.boost` diagnostic, which
+is documented as a stochastic sampling gap (R RNG-driven CV-index, feature
+subset, and up/down-sampling draws cannot be reproduced bit-for-bit with
+NumPy's RNG). No deterministic parity gap is excluded from the suite.
+
+### Manual comparison harness
+
+No `scripts/compare_nns.py` / `scripts/compare_nns_r13.py` diagnostic scripts
+exist in this repository, so there was no `$RPM`-based univariate `NNS.reg`
+extraction to fix. The univariate regression-point diagnostic lives in
+`tests/parity/test_r13_smoke.py` and already extracts regression points via
+`NNS.reg(..., multivariate.call = TRUE)$y`, matching the supported R-side
+extraction.
+
+### ARMA nonseasonal nonlinear reconciliation
+
+The suspected live mismatch (Python `[125.25, 107.75, 158.75, 213.66]` vs R
+`[128.50, 113.50, 155.50, 213.66]`) was re-run against live vendored R NNS
+13.0 before regeneration:
+
+- Live R `NNS.ARMA(series, h = 4, seasonal.factor = FALSE, method = "nonlin")`
+  on the 24-point AirPassengers-style series returned
+  `128.5, 113.5, 155.5, 213.6666666667`.
+- Python `nns_arma(series, h=4, seasonal_factor=False, method="nonlin")`
+  returned `128.5, 113.5, 155.5, 213.66666667`.
+
+Python matches live R 13.0 exactly; the previously reported divergent Python
+values could not be reproduced with the current implementation. The hardcoded
+smoke expectations in `tests/parity/test_r13_smoke.py` (LPM/UPM, copulas,
+regression points, seasonal ARMA, nonseasonal nonlinear ARMA, seeded stack)
+were each verified against this fresh live R 13.0 run and required no changes.

scripts/regenerate_r_cache.py

@@ -3,21 +3,43 @@
 
 CI should not run this script. It intentionally clears cache-only/offline toggles
 and invokes pytest so tests/_r.py can refresh tests/_r_cache.json as needed.
+
+Usage::
+
+    python scripts/regenerate_r_cache.py [--fresh] [-- PYTEST_ARGS...]
+
+By default, existing cache entries are reused and only cache misses call R.
+With ``--fresh``, the existing ``tests/_r_cache.json`` is moved aside to
+``tests/_r_cache.json.bak`` and regeneration starts from an empty cache, so
+every entry is produced by a live R call. ``--fresh`` refuses to run in CI and
+requires a working local R NNS install reporting ``packageVersion("NNS")`` of
+13.0 (see ``scripts/install_local_r_nns.py``).
 """
 
 from __future__ import annotations
 
+import argparse
 import json
 import os
+import shutil
 import subprocess
 import sys
 from pathlib import Path
 from typing import Any
 
 _CACHE_PATH = Path(__file__).resolve().parents[1] / "tests" / "_r_cache.json"
+_CACHE_BACKUP_PATH = _CACHE_PATH.with_suffix(".json.bak")
 _NNS_VERSION = "13.0"
 _SCHEMA_VERSION = 1
 
+_OFFLINE_TOGGLES = (
+    "PYNNS_R_CACHE_ONLY",
+    "NNS_R_CACHE_ONLY",
+    "PYNNS_OFFLINE",
+    "NNS_OFFLINE",
+    "CI",
+)
+
 
 def _validate_cache() -> int:
     if not _CACHE_PATH.exists():
@@ -71,15 +93,95 @@ def _validate_cache() -> int:
         )
         return 1
 
+    print(f"OK: {_CACHE_PATH} contains {len(entries)} entries for NNS {_NNS_VERSION}.")
+    return 0
+
+
+def _running_in_ci() -> bool:
+    return os.environ.get("CI", "").lower() in {"1", "true", "yes"} or bool(
+        os.environ.get("GITHUB_ACTIONS")
+    )
+
+
+def _verify_live_r_nns() -> int:
+    """Confirm a local R NNS install reports the expected version before a fresh run."""
+
+    rscript = shutil.which("Rscript")
+    if rscript is None:
+        print(
+            "ERROR: --fresh requires Rscript on PATH; "
+            "run scripts/install_local_r_nns.py first.",
+            file=sys.stderr,
+        )
+        return 1
+    probe = subprocess.run(
+        [
+            rscript,
+            "-e",
+            "suppressPackageStartupMessages(library(NNS)); "
+            "cat(as.character(packageVersion('NNS')))",
+        ],
+        capture_output=True,
+        text=True,
+        check=False,
+    )
+    if probe.returncode != 0:
+        print(
+            "ERROR: --fresh could not load R NNS:\n" + probe.stderr,
+            file=sys.stderr,
+        )
+        return 1
+    version = probe.stdout.strip()
+    if version != _NNS_VERSION:
+        print(
+            f"ERROR: --fresh requires R NNS {_NNS_VERSION!r}; "
+            f"installed version is {version!r}.",
+            file=sys.stderr,
+        )
+        return 1
+    print(f"OK: live R NNS {version} detected for fresh regeneration.")
     return 0
 
 
 def main() -> int:
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument(
+        "--fresh",
+        action="store_true",
+        help=(
+            "Move the existing cache to tests/_r_cache.json.bak and regenerate every "
+            "entry from a live R call. Refuses to run in CI."
+        ),
+    )
+    parser.add_argument(
+        "pytest_args",
+        nargs="*",
+        help="Arguments passed to pytest after an optional '--' separator.",
+    )
+    parsed = parser.parse_args()
+
+    if parsed.fresh:
+        if _running_in_ci():
+            print(
+                "ERROR: --fresh must not run in CI; it deletes the committed cache "
+                "and requires a local R NNS install.",
+                file=sys.stderr,
+            )
+            return 1
+        verify_status = _verify_live_r_nns()
+        if verify_status:
+            return verify_status
+        if _CACHE_PATH.exists():
+            _CACHE_PATH.replace(_CACHE_BACKUP_PATH)
+            print(f"Moved existing cache to {_CACHE_BACKUP_PATH}; starting from empty cache.")
+        else:
+            print("No existing cache found; starting from empty cache.")
+
     env = os.environ.copy()
-    for name in ("PYNNS_R_CACHE_ONLY", "PYNNS_OFFLINE", "CI"):
+    for name in _OFFLINE_TOGGLES:
         env.pop(name, None)
 
-    args = sys.argv[1:]
+    args = parsed.pytest_args
     if args[:1] == ["--"]:
         args = args[1:]
     if not args:

src/nns/arma.py

@@ -753,8 +753,10 @@ def _numeric_seasonal_weights(
     lags: NDArray[np.int64],
 ) -> NDArray[np.float64]:
     output = np.empty(lags.size, dtype=np.float64)
-    for index, lag in enumerate(lags):
-        rev_var = variable[:: -int(lag)]
+    for index in range(lags.size):
+        # R NNS 13.0 reverses by the candidate's POSITION in seasonal.factor
+        # (variable[seq(length(variable), 1, -i)]), not by its lag value.
+        rev_var = variable[:: -(index + 1)]
         with np.errstate(invalid="ignore", divide="ignore"):
             output[index] = abs(
                 np.float64(np.std(rev_var, ddof=1)) / np.float64(np.mean(rev_var))

src/nns/boost.py

@@ -282,9 +282,11 @@ def _nns_boost_core(
         type=type_value,
         pred_int=pred_int,
     )
-    estimates = np.asarray(final_fit["stack"], dtype=np.float64)
-    if estimates.size == 0 or np.any(np.isnan(estimates)):
+    stack_estimates = final_fit["stack"]
+    if stack_estimates is None or np.asarray(stack_estimates, dtype=np.float64).size == 0:
         estimates = np.asarray(final_fit["reg"], dtype=np.float64)
+    else:
+        estimates = np.asarray(stack_estimates, dtype=np.float64)
     estimates = _fill_nan_with_gravity(estimates)
     if type_value == "class":
         estimates = _round_clamp_classes(estimates, y_train)
@@ -294,7 +296,6 @@ def _nns_boost_core(
         "pred.int": final_fit["pred.int"],
         "feature.weights": weights[order_idx],
         "feature.frequency": counts[order_idx],
-        "n.best": final_fit["NNS.reg.n.best"],
     }