Wire plot=True on compute functions to render via nns.plotting

Previously the compute functions accepted R-compatible plot arguments and
immediately `del`'d them, so `plot=True` produced nothing — a strict validator
expecting a figure would fail. Now `plot=True` renders a Matplotlib figure as a
side effect through the nns.plotting layer while the value-only return contract
(which the parity suite depends on) is unchanged.

Wired:
- nns_reg: plot / plot_regions -> plot_nns_reg; residual_plot -> residual scatter
- nns_m_reg: plot -> fitted-vs-actual; residual_plot -> residual scatter
- nns_arma / nns_arma_optim -> plot_nns_arma / plot_nns_arma_optim
- nns_cdf (univariate) -> plot_nns_cdf
- nns_seas -> plot_nns_seas (compute split into _nns_seas_compute)

Rendering only fires for results that carry the needed series, and default
plot=False opens no figure (verified). Adds tests/plotting/test_compute_plot_flag.py
asserting plot=True creates a figure and returns the same value as plot=False.
Updates the plot-parity policy and README.

Author: claude

README.md

@@ -142,7 +142,7 @@ Important boundaries:
 - Stochastic exact stream parity is not expected because Python paths use NumPy random generation.
 - Factor and class ordering should be passed explicitly when ordering matters.
 - Direct raw-factor `nns_m_reg(..., factor_2_dummy=True)` is intentionally guarded. Use `prepare_factor_predictors(...)` before `nns_m_reg(...)`.
-- Compute functions' `plot` arguments are ignored and data is returned instead; visual plotting is a separate API in `nns.plotting`, color/element-faithful to R but not pixel-diffed.
+- Compute functions still return values, not figures; passing `plot=True` (where R has it) additionally renders a Matplotlib figure as a side effect via the `nns.plotting` layer, which is color/element-faithful to R but not pixel-diffed. The plot functions can also be called directly on a computed result.
 
 See [behavior conventions](docs/conventions.md) for detailed compatibility notes.
 

docs/plot_parity_policy.md

@@ -50,9 +50,13 @@ colors and which element they sit on*, never rendered images.
   a plotting API, so Python parity calls pass the R `plot = FALSE` equivalent
   and assert only on returned values.
 
-When a ported function has an R `plot` argument, the Python API either omits the
-argument entirely or treats plotting as out of scope; only the value-bearing
-return is asserted in parity tests.
+When a ported function has an R `plot` argument, the Python function keeps its
+value-only **return** contract (parity asserts only on the returned value). As a
+side effect, passing `plot=True` renders a Matplotlib figure through the
+`nns.plotting` layer — `nns_reg`, `nns_m_reg`, `nns_arma`, `nns_arma_optim`,
+`nns_cdf`, and `nns_seas` are wired this way (plus `residual_plot=True` for the
+regression functions). The figures are color/element-faithful but never
+pixel-compared, and computation with the default `plot=False` opens no figure.
 
 ## Inventory of committed graphics artifacts
 

src/nns/arma.py

@@ -32,7 +32,7 @@ def nns_arma_optim(
     plot: bool = False,
 ) -> dict[str, Any]:
     """Optimize seasonal factors for :func:`nns_arma` like R's ``NNS.ARMA.optim``."""
-    del ncores, print_trace, plot
+    del ncores, print_trace
 
     values = _as_variable(variable)
     original_values = values.copy()
@@ -326,7 +326,7 @@ def nns_arma_optim(
         lower_pi = np.maximum(0.0, lower_pi)
         upper_pi = np.maximum(0.0, upper_pi)
 
-    return {
+    result = {
         "periods": nns_periods,
         "weights": nns_weights,
         "obj.fn": nns_score,
@@ -339,6 +339,11 @@ def nns_arma_optim(
         "lower.pred.int": lower_pi,
         "upper.pred.int": upper_pi,
     }
+    if plot:
+        from nns.plotting.arma import plot_nns_arma_optim
+
+        plot_nns_arma_optim(result, original_values)
+    return result
 
 
 def nns_arma(
@@ -360,12 +365,22 @@ def nns_arma(
     random_seed: int | None = None,
 ) -> NDArray[np.float64] | dict[str, NDArray[np.float64]]:
     """Autoregressive NNS forecast matching R's installed NNS.ARMA behavior."""
-    del plot, seasonal_plot
+    del seasonal_plot
 
     horizon = int(h)
     if horizon < 1:
         raise ValueError("h must be a positive integer.")
     values = _as_variable(variable)
+
+    def _finish(
+        forecast: NDArray[np.float64] | dict[str, NDArray[np.float64]],
+    ) -> NDArray[np.float64] | dict[str, NDArray[np.float64]]:
+        if plot:
+            from nns.plotting.arma import plot_nns_arma
+
+            ts = int(training_set) if training_set is not None else int(values.size)
+            plot_nns_arma(forecast, values, training_set=ts)
+        return forecast
     if _is_numeric_seasonal(seasonal_factor) and dynamic:
         raise ValueError(
             'Hmmm...Seems you have "seasonal.factor" specified and "dynamic = TRUE".  '
@@ -389,12 +404,12 @@ def nns_arma(
 
     estimates = np.zeros(horizon, dtype=np.float64)
     if not _is_numeric_seasonal(seasonal_factor) and np.ptp(values) == 0.0:
-        return _with_prediction_intervals(
+        return _finish(_with_prediction_intervals(
             estimates,
             lin_residual=0.0,
             pred_int=pred_int,
             random_seed=random_seed,
-        )
+        ))
     lags, lag_weights = _resolve_lags_and_weights(
         values,
         seasonal_factor=seasonal_factor,
@@ -416,7 +431,7 @@ def nns_arma(
             method=method_l,
             shrink=shrink,
         )
-        return estimates
+        return _finish(estimates)
 
     current = values
     lin_regression_estimates = np.array([], dtype=np.float64)
@@ -490,12 +505,12 @@ def nns_arma(
         if not np.isfinite(lin_resid):
             lin_resid = 0.0
 
-    return _with_prediction_intervals(
+    return _finish(_with_prediction_intervals(
         estimates,
         lin_residual=lin_resid,
         pred_int=pred_int,
         random_seed=random_seed,
-    )
+    ))
 
 
 def _valid_arma_optim_seasonals(

src/nns/cdf.py

@@ -20,7 +20,6 @@ def nns_cdf(
     names: Sequence[str] | None = None,
 ) -> dict[str, object]:
     """Partial-moment CDF wrapper matching R's non-plotting NNS.CDF paths."""
-    del plot
     type_value = type.lower()
     if type_value not in {"cdf", "survival", "hazard", "cumulative hazard"}:
         raise ValueError("invalid type")
@@ -29,12 +28,28 @@ def nns_cdf(
     if values.ndim == 0:
         values = values.reshape(1)
     if values.ndim == 1 or (values.ndim == 2 and values.shape[1] == 1):
-        return _univariate_cdf(values.reshape(-1), float(degree), target, type_value)
+        result = _univariate_cdf(values.reshape(-1), float(degree), target, type_value)
+        if plot:
+            _render_cdf(result, target)
+        return result
     if values.ndim == 2:
+        # Multivariate CDF has no faithful single-Axes plot; plot is a no-op here.
         return _multivariate_cdf(values, float(degree), target, type_value, names)
     raise ValueError("variable must be a vector or 2D matrix.")
 
 
+def _render_cdf(result: dict[str, object], target: float | NDArray[np.float64] | None) -> None:
+    """Render the univariate NNS.CDF figure as a side effect of ``plot=True``."""
+    from nns.plotting.partial_moments import plot_nns_cdf
+
+    plot_target: float | None = None
+    if target is not None:
+        coords = np.asarray(target, dtype=np.float64).reshape(-1)
+        if coords.size == 1:
+            plot_target = float(coords[0])
+    plot_nns_cdf(result, target=plot_target)
+
+
 def _univariate_cdf(
     values: NDArray[np.float64],
     degree: float,

src/nns/multivariate_regression.py

@@ -39,8 +39,8 @@ def nns_m_reg(
     confidence_interval: float | None = None,
     class_levels: list[object] | None = None,
 ) -> MRegResult:
-    """Multivariate numeric regression matching R's non-plotting NNS.M.reg path."""
-    del plot, residual_plot, location, dist, return_values, plot_regions, ncores
+    """Multivariate numeric regression matching R's NNS.M.reg path."""
+    del location, dist, return_values, plot_regions, ncores
     type_value = _normalize_type(type)
     x_values, y_values = _validate_inputs(
         x,
@@ -107,14 +107,50 @@ def nns_m_reg(
         confidence_interval=confidence_interval,
     )
     r2 = _class_accuracy(y_values, fitted_y) if type_value == "class" else _r2(y_values, fitted_y)
-    return {
+    result: MRegResult = {
         "R2": r2,
         "rhs.partitions": _rhs_partitions_dict(reg_points_matrix),
         "RPM": _rpm_dict(rpm),
         "Point.est": _point_output(point_predictions),
         "pred.int": pred_int,
         "Fitted.xy": fitted,
     }
+    if plot or residual_plot:
+        _render_m_reg(fitted, plot=plot, residual_plot=residual_plot)
+    return result
+
+
+def _render_m_reg(
+    fitted: dict[str, NDArray[np.float64] | NDArray[np.str_]],
+    *,
+    plot: bool,
+    residual_plot: bool,
+) -> None:
+    """Render multivariate fitted-vs-actual / residual diagnostics on ``plot=True``.
+
+    The synthetic predictors make a single x-axis ill-defined, so ``plot`` shows
+    fitted vs actual (steelblue points, red 1:1 line) and ``residual_plot`` shows
+    residuals about zero -- a figure is still produced, faithful to R's colors.
+    """
+    from nns.plotting._mpl import resolve_ax
+
+    y = np.asarray(fitted["y"], dtype=np.float64)
+    y_hat = np.asarray(fitted["y.hat"], dtype=np.float64)
+    if plot and y.size and y.size == y_hat.size:
+        ax = resolve_ax(None)
+        ax.scatter(y, y_hat, color="steelblue")
+        lo, hi = float(min(y.min(), y_hat.min())), float(max(y.max(), y_hat.max()))
+        ax.plot([lo, hi], [lo, hi], color="red")
+        ax.set_xlabel("y")
+        ax.set_ylabel("y.hat")
+        ax.set_title("NNS.M.reg Fitted vs Actual")
+    if residual_plot:
+        residuals = np.asarray(fitted.get("residuals", []), dtype=np.float64)
+        if residuals.size:
+            ax = resolve_ax(None)
+            ax.scatter(np.arange(1, residuals.size + 1), residuals, color="steelblue")
+            ax.axhline(0.0, color="red")
+            ax.set_title("NNS.M.reg Residual Plot")
 
 
 def _validate_inputs(

src/nns/regression.py

@@ -57,10 +57,10 @@ def nns_reg(
     factor_levels: Sequence[object] | Sequence[Sequence[object] | None] | None = None,
 ) -> dict[str, Any]:
     """Univariate numeric port of R's NNS.reg."""
-    del return_values, plot, plot_regions, residual_plot, ncores
+    del return_values, ncores
 
     if dim_red_method is not None:
-        return _nns_reg_dimred(
+        result = _nns_reg_dimred(
             x,
             y,
             factor_2_dummy=factor_2_dummy,
@@ -80,6 +80,11 @@ def nns_reg(
             class_levels=class_levels,
             factor_levels=factor_levels,
         )
+        _maybe_render_reg(
+            result, plot=plot, plot_regions=plot_regions,
+            residual_plot=residual_plot, point_est=point_est,
+        )
+        return result
 
     type_value = _normalize_type(type)
     if type_value == "class":
@@ -101,7 +106,7 @@ def nns_reg(
         dispatch_n_best = n_best
         if type_value == "class" and dispatch_n_best is None:
             dispatch_n_best = 1
-        return nns_m_reg(
+        result = nns_m_reg(
             np.asarray(x_for_dispatch, dtype=np.float64),
             y_matrix_values,
             factor_2_dummy=False,
@@ -117,6 +122,11 @@ def nns_reg(
             confidence_interval=confidence_interval,
             class_levels=class_levels,
         )
+        _maybe_render_reg(
+            result, plot=plot, plot_regions=plot_regions,
+            residual_plot=residual_plot, point_est=point_est,
+        )
+        return result
 
     del tau, threshold, n_best, dist
     x_values, y_values = _validate_univariate_inputs(
@@ -137,7 +147,7 @@ def nns_reg(
     )
     noise = _validate_noise_reduction(noise_reduction)
     point_values = _as_point_est(point_for_dispatch)
-    return _nns_reg_univariate_core(
+    result = _nns_reg_univariate_core(
         x_values,
         y_values,
         order=order,
@@ -150,6 +160,47 @@ def nns_reg(
         equation=None,
         x_star=None,
     )
+    _maybe_render_reg(
+        result, plot=plot, plot_regions=plot_regions,
+        residual_plot=residual_plot, point_est=point_est,
+    )
+    return result
+
+
+def _maybe_render_reg(
+    result: dict[str, Any],
+    *,
+    plot: bool,
+    plot_regions: bool,
+    residual_plot: bool,
+    point_est: NDArray[np.float64] | float | None,
+) -> None:
+    """Render the NNS.reg figure(s) as a side effect when a plot flag is set.
+
+    Plotting is decoupled from computation: this only fires for the standard
+    univariate result (one that carries ``Fitted.xy``) so the value-only return
+    contract is unchanged. ``plot``/``plot_regions`` draw the regression figure;
+    ``residual_plot`` draws a residual scatter.
+    """
+    if not (plot or plot_regions or residual_plot):
+        return
+    fitted = result.get("Fitted.xy") if isinstance(result, dict) else None
+    if not isinstance(fitted, dict) or "x" not in fitted:
+        return
+    if plot or plot_regions:
+        from nns.plotting.regression import plot_nns_reg
+
+        plot_nns_reg(result, point_est=point_est)
+    if residual_plot:
+        from nns.plotting._mpl import resolve_ax
+
+        xs = np.asarray(fitted["x"], dtype=np.float64)
+        residuals = np.asarray(fitted.get("residuals", []), dtype=np.float64)
+        if residuals.size and residuals.size == xs.size:
+            ax = resolve_ax(None)
+            ax.scatter(xs, residuals, color="steelblue")
+            ax.axhline(0.0, color="red")
+            ax.set_title("NNS Residual Plot")
 
 
 def prepare_factor_predictors(

src/nns/seasonality.py

@@ -19,8 +19,21 @@ def nns_seas(
     mod_only: bool = True,
     plot: bool = False,
 ) -> SeasonalityResult:
-    """Seasonality test matching R's NNS.seas non-plotting path."""
-    del plot
+    """Seasonality test matching R's NNS.seas; ``plot=True`` renders the figure."""
+    result = _nns_seas_compute(variable, modulo=modulo, mod_only=mod_only)
+    if plot:
+        from nns.plotting.seasonality import plot_nns_seas
+
+        plot_nns_seas(result)
+    return result
+
+
+def _nns_seas_compute(
+    variable: NDArray[np.float64],
+    *,
+    modulo: int | list[int] | NDArray[np.int64] | None = None,
+    mod_only: bool = True,
+) -> SeasonalityResult:
     values = _validate_variable(variable)
     modulo_values = None if modulo is None else _as_modulo(modulo)
     cache_key = _cache_key(values, modulo_values, mod_only)

tests/plotting/conftest.py

@@ -5,3 +5,7 @@
 import matplotlib
 
 matplotlib.use("Agg")
+
+# Plotting tests open many short-lived figures (closed per-test); don't warn
+# about the open-figure count if this session shares a process with others.
+matplotlib.rcParams["figure.max_open_warning"] = 0