feat(gam): close general-family multismooth surfaces

nampy/gam/diagnostics/concurvity.py

@@ -15,6 +15,7 @@
     _term_full_coefficient_indices,
 )
 from ..predict.linear_predictor_matrix import build_lpmatrix
+from ..term_labels import mgcv_term_display_label
 
 
 def _term_indices_for_concurvity(model, n_coef: int):
@@ -37,7 +38,7 @@ def _term_indices_for_concurvity(model, n_coef: int):
         if idx.size == 0:
             continue
         smooth_starts.append(int(np.min(idx)))
-        blocks.append((str(tb.label), idx))
+        blocks.append((mgcv_term_display_label(tb), idx))
 
     if len(blocks) == 0:
         raise ValueError("No smooth or parametric components available for concurvity.")

nampy/gam/diagnostics/derivatives.py

@@ -12,6 +12,7 @@
     _term_blocks_seq,
     _term_full_coefficient_indices,
 )
+from ..term_labels import mgcv_term_display_label
 
 
 @dataclass(frozen=True)
@@ -65,7 +66,7 @@ def smooth_derivative(model, *, X=None, smooth_number: int = 1, deriv: int = 1):
         derivative=values,
         se=np.sqrt(np.maximum(variance, 0.0)),
         derivative_matrix=Xd,
-        term_label=str(term.label),
+        term_label=mgcv_term_display_label(term),
         order=deriv,
     )
 

nampy/gam/diagnostics/k_check.py

@@ -13,6 +13,7 @@
     _term_blocks_seq,
 )
 from ..predict.linear_predictor_matrix import build_lpmatrix
+from ..term_labels import mgcv_term_display_label
 from .residuals import residuals_gam
 
 
@@ -177,7 +178,7 @@ def k_check(model, subsample: int = 5000, n_rep: int = 400, seed: int | None = N
         edf_by_term = edf_all
     perm_col = 0
     for i, tb in enumerate(term_blocks):
-        label = str(tb.label)
+        label = mgcv_term_display_label(tb)
         X_term = _numeric_feature_block(model, tb, row_idx)
         k_prime = int(tb.coef_slice.stop - tb.coef_slice.start)
         edf = float(edf_by_term[i])

nampy/gam/diagnostics/plots.py

@@ -36,7 +36,7 @@
 from ..predict.predictions import _term_has_absorbed_constraint
 from ..predict.terms import _prediction_term_groups
 from ..smooths.categorical import factor_levels_from_metadata
-from ..term_labels import normalize_mgcv_term_label
+from ..term_labels import mgcv_term_display_label
 from .residuals import _prior_weights
 
 __all__ = [
@@ -462,7 +462,8 @@ def prepare_plot_gam_data(
 
     pd_list = []
     for i, tb in enumerate(smooth_blocks):
-        label = _sub_edf(str(tb.label), edf_map.get(id(tb), float("nan")))
+        display_label = mgcv_term_display_label(tb)
+        label = _sub_edf(display_label, edf_map.get(id(tb), float("nan")))
         P = _prepare_smooth(
             model,
             tb,
@@ -535,7 +536,7 @@ def prepare_plot_gam_data(
             P["se"] = False
 
         if partial_resids:
-            normalized = str(normalize_mgcv_term_label(str(tb.label)))
+            normalized = display_label
             group_index = None
             for gi, glabel in enumerate(term_group_labels):
                 if glabel == normalized:

nampy/gam/fit/selection/postfit.py

@@ -12,7 +12,9 @@
     _n_smoothing_params,
     _penalty_blocks_seq,
     _require_fitted,
+    _term_blocks_seq,
 )
+from ...term_labels import multi_predictor_term_label
 from .criteria.dispatch import criterion_gradient, criterion_hessian, criterion_value
 from .criteria.ml_reml import resolve_ml_reml_scoring_backend
 
@@ -131,6 +133,7 @@ def _gam_vcomp_names(model) -> list[str]:
         return []
 
     names: list[str | None] = [None] * n_sp
+    terms = list(_term_blocks_seq(model))
     for pb in _penalty_blocks_seq(model):
         idx = int(getattr(pb, "smoothing_index", -1))
         if idx < 0 or idx >= n_sp or names[idx] is not None:
@@ -141,7 +144,16 @@ def _gam_vcomp_names(model) -> list[str]:
             label = meta.get("label", None)
         if label is None:
             label = getattr(pb, "label", None)
-        names[idx] = _normalize_vcomp_label(label)
+        normalized = _normalize_vcomp_label(label)
+        term_index = int(getattr(pb, "term_index", -1))
+        if normalized is not None and 0 <= term_index < len(terms):
+            owner = terms[term_index]
+            normalized = multi_predictor_term_label(
+                normalized,
+                predictor_index=int(getattr(owner, "predictor_index", 0)),
+                term_type=str(getattr(owner, "term_type", "")),
+            )
+        names[idx] = normalized
 
     return [name if name is not None else f"sp_{i}" for i, name in enumerate(names)]
 

nampy/gam/inference/anova.py

@@ -34,6 +34,7 @@
     _term_blocks_seq,
     _term_full_coefficient_indices,
 )
+from ..term_labels import mgcv_term_display_label
 from .chi_square_mixtures import psum_chisq
 
 
@@ -45,9 +46,7 @@ def _scale_estimated(model) -> bool:
 
 
 def _formula_term_label(tb) -> str:
-    metadata = dict(getattr(tb, "metadata", {}) or {})
-    formula_term = metadata.get("formula_term", None)
-    return str(getattr(tb, "label", "")) if formula_term is None else str(formula_term)
+    return mgcv_term_display_label(tb, formula_parametric=True)
 
 
 def _parametric_term_groups(model):
@@ -56,7 +55,7 @@ def _parametric_term_groups(model):
         if str(getattr(tb, "term_type", "")) != "parametric":
             continue
         label = _formula_term_label(tb)
-        key = ("parametric", label)
+        key = ("parametric", int(getattr(tb, "predictor_index", 0)), label)
         if groups and groups[-1]["key"] == key:
             groups[-1]["blocks"].append(tb)
             continue
@@ -620,7 +619,7 @@ def _term_table(
         )
         smooth_rows.append(
             {
-                "label": str(tb.label),
+                "label": mgcv_term_display_label(tb),
                 "edf": edf_i,
                 "ref_df": ref_df,
                 "wald_stat": stat_out,

nampy/gam/inference/summary.py

@@ -32,6 +32,7 @@
     _term_blocks_seq,
     _term_full_coefficient_indices,
 )
+from ..term_labels import mgcv_term_display_label
 from .anova import (
     _parametric_term_groups,
     _residual_df,
@@ -116,12 +117,13 @@ def _parametric_coefficient_indices(model) -> tuple[list[int], list[str]]:
     for group in _parametric_term_groups(model):
         for tb in group["blocks"]:
             full_indices = _term_full_coefficient_indices(model, tb)
+            term_label = mgcv_term_display_label(tb)
             for j, full_index in enumerate(full_indices):
                 indices.append(int(full_index))
                 names.append(
-                    str(tb.label)
+                    term_label
                     if full_indices.size == 1
-                    else f"{tb.label}.{j}"
+                    else f"{term_label}.{j}"
                 )
     return indices, names
 

nampy/gam/predict/terms.py

@@ -17,7 +17,7 @@
     _term_blocks_seq,
     _term_full_coefficient_indices,
 )
-from ..term_labels import normalize_mgcv_term_label
+from ..term_labels import multi_predictor_term_label, normalize_mgcv_term_label
 
 
 def _parametric_formula_term(term) -> str | None:
@@ -26,18 +26,6 @@ def _parametric_formula_term(term) -> str | None:
     return None if formula_term is None else str(formula_term)
 
 
-def _multi_predictor_term_label(label: str, *, predictor_index: int, term_type: str):
-    """Apply mgcv's formula-list suffix to later-predictor term labels."""
-    if predictor_index <= 0:
-        return label
-    if term_type == "parametric":
-        return f"{label}.{predictor_index}"
-    open_index = label.find("(")
-    if open_index < 0:
-        return f"{label}.{predictor_index}"
-    return f"{label[:open_index]}.{predictor_index}{label[open_index:]}"
-
-
 def _prediction_term_groups(model):
     """Return ordered mgcv term groups without treating labels as identity."""
     groups = []
@@ -47,7 +35,7 @@ def _prediction_term_groups(model):
         predictor_name = str(getattr(term, "predictor_name", "predictor_0"))
         if term_type == "parametric":
             formula_term = _parametric_formula_term(term)
-            group_label = _multi_predictor_term_label(
+            group_label = multi_predictor_term_label(
                 formula_term or str(getattr(term, "label", "")),
                 predictor_index=predictor_index,
                 term_type=term_type,
@@ -76,7 +64,7 @@ def _prediction_term_groups(model):
             )
             continue
 
-        group_label = _multi_predictor_term_label(
+        group_label = multi_predictor_term_label(
             str(normalize_mgcv_term_label(getattr(term, "label", ""))),
             predictor_index=predictor_index,
             term_type=term_type,

nampy/gam/results/snapshots.py

@@ -37,7 +37,7 @@
     _term_full_coefficient_indices,
 )
 from ..predict.terms import _prediction_term_groups
-from ..term_labels import normalize_mgcv_term_label
+from ..term_labels import mgcv_term_display_label, normalize_mgcv_term_label
 
 
 def _as_pred_or_scalar_array(value):
@@ -430,7 +430,7 @@ def build_parity_snapshot(model, X=None, include_covariances=False):
                 continue
             smooth_blocks.append(tb)
             full_idx = _term_full_coefficient_indices(core, tb)
-            smooth_labels.append(_normalize_reference_term_label(tb.label))
+            smooth_labels.append(mgcv_term_display_label(tb))
             if _cov_bayes(core) is not None:
                 smooth_cov_bayes.append(
                     np.asarray(

nampy/gam/term_labels.py

@@ -84,4 +84,38 @@ def normalize_mgcv_term_label(label):
     return text
 
 
-__all__ = ["normalize_mgcv_term_label"]
+def multi_predictor_term_label(
+    label: str, *, predictor_index: int, term_type: str
+) -> str:
+    """Apply mgcv's formula-list suffix to later-predictor term labels."""
+    if int(predictor_index) <= 0:
+        return str(label)
+    if str(term_type) == "parametric":
+        return f"{label}.{int(predictor_index)}"
+    open_index = str(label).find("(")
+    if open_index < 0:
+        return f"{label}.{int(predictor_index)}"
+    return f"{label[:open_index]}.{int(predictor_index)}{label[open_index:]}"
+
+
+def mgcv_term_display_label(term, *, formula_parametric: bool = False) -> str:
+    """Return a compiled term's predictor-aware public mgcv label."""
+    term_type = str(getattr(term, "term_type", ""))
+    label = str(getattr(term, "label", ""))
+    if term_type == "parametric" and formula_parametric:
+        metadata = dict(getattr(term, "metadata", {}) or {})
+        label = str(metadata.get("formula_term", label))
+    elif term_type != "parametric":
+        label = str(normalize_mgcv_term_label(label))
+    return multi_predictor_term_label(
+        label,
+        predictor_index=int(getattr(term, "predictor_index", 0)),
+        term_type=term_type,
+    )
+
+
+__all__ = [
+    "mgcv_term_display_label",
+    "multi_predictor_term_label",
+    "normalize_mgcv_term_label",
+]

tests/SUBSYSTEM_COVERAGE.md

@@ -25,7 +25,7 @@ local development notes rather than duplicated here.
 | ML/REML backend routing | `nampy/gam/fit/selection/criteria/ml_reml.py` | `tests/optimization/test_gam_owner_routing_objective_contracts.py`, `tests/optimization/test_mgcv_gaussian_backend_selection.py` | Exact vs dynamic vs PIRLS vs general-family selection. |
 | Objective wrappers / optimizer wiring | `nampy/gam/fit/selection/optimize/objectives.py`, `.../driver.py` | `tests/optimization/test_gam_owner_routing_objective_contracts.py`, `tests/optimization/test_mgcv_parametric_only_parity.py`, `tests/optimization/test_mgcv_outer_optimization_parity.py`, `tests/optimization/test_mgcv_optimization_lifecycle_parity.py` | Owner contracts first, direct empty-smoothing-vector parity, mgcv trace parity, then lifecycle parity. |
 | Postfit smoothing diagnostics | `nampy/gam/fit/selection/postfit.py` | `tests/optimization/test_gam_postfit_owner_contracts.py`, `tests/optimization/test_mgcv_vcomp_parity.py`, `tests/optimization/test_mgcv_sp_vcov_stage_parity.py` | Endpoint diagnostics, Hessian sourcing, smoothing covariance surfaces, and stage-local `sp.vcov` / unconditional-covariance checkpoints. |
-| General-family fixed-smoothing / postprocess | `nampy/gam/fit/solvers/general_family/fixed_smoothing.py`, `.../newton.py` | `tests/families/test_gam_general_family_owner_contracts.py`, `tests/optimization/test_mgcv_fixed_inner_fit_parity.py`, `tests/optimization/test_mgcv_general_family_preoptimization_parity.py`, `tests/families/test_general_family_mgcv_parity.py`, `tests/parity/test_mgcv_under_tested_supported_combinations.py` | Owner precedence plus mgcv `gam.fit5` parity; reparameterized and original-coordinate singleton/multi-penalty `Sl` blocks are covered through setup, derivatives, roots/totals, and full fits. Structured `re`, `fs`, and linked-`sz` cases include fixed and optimized `fs` behavior and an `fs` block in linear predictor two. |
+| General-family fixed-smoothing / postprocess | `nampy/gam/fit/solvers/general_family/fixed_smoothing.py`, `.../newton.py` | `tests/families/test_gam_general_family_owner_contracts.py`, `tests/optimization/test_mgcv_fixed_inner_fit_parity.py`, `tests/optimization/test_mgcv_general_family_preoptimization_parity.py`, `tests/families/test_general_family_mgcv_parity.py`, `tests/parity/test_mgcv_snapshot_extended_matrix.py`, `tests/parity/test_mgcv_under_tested_supported_combinations.py` | Owner precedence plus mgcv `gam.fit5` parity; reparameterized and original-coordinate singleton/multi-penalty `Sl` blocks are covered through setup, derivatives, roots/totals, and full fits. A four-smooth/two-predictor gaulss case covers wrapped blocks, `sp.vcov`, inference, and diagnostics. Structured `re`, `fs`, and linked-`sz` cases include fixed and optimized `fs` behavior and an `fs` block in linear predictor two. |
 | Diagnostics owners | `nampy/gam/diagnostics/residuals.py`, `concurvity.py`, `summary.py`, `plots.py` | `tests/diagnostics/test_gam_diagnostics_owner_contracts.py`, `tests/diagnostics/test_gam_plot_and_public_api_contracts.py`, `tests/parity/test_mgcv_secondary_diagnostics_parity.py` | Owner-level residual/summary/plot contracts plus direct secondary-diagnostics parity. |
 | Prediction / inference / diagnostics | `nampy/gam/predict/`, `nampy/gam/inference/`, `nampy/gam/diagnostics/` | `tests/parity/test_mgcv_output_parity.py`, `tests/parity/test_mgcv_prediction_arguments_parity.py`, `tests/parity/test_mgcv_prediction_inference_diagnostics_parity.py`, `tests/parity/test_mgcv_general_family_lpmatrix_stage_parity.py`, `tests/parity/test_mgcv_general_family_prediction_stage_parity.py`, `tests/parity/test_mgcv_inference_stage_parity.py`, `tests/diagnostics/test_mgcv_general_family_secondary_diagnostics_parity.py` | Public-surface parity plus direct `block.size`, `newdata.guaranteed`, `na.action`, `unconditional`, `iterms.type`, and stage-local general-family checkpoints. |
 | Parity snapshot / trace tooling | `nampy/gam/parity/` | `tests/parity/test_gam_parity_owner_contracts.py`, `tests/parity/test_gam_results_api_stage_owner_contracts.py`, `tests/optimization/test_mgcv_score_hist_trace_parity.py`, `tests/optimization/test_mgcv_outer_optimization_parity.py`, `tests/optimization/test_mgcv_optimization_lifecycle_parity.py`, `tests/optimization/test_mgcv_inner_trace_parity.py`, `tests/optimization/test_mgcv_joint_branch_trace_parity.py` | Localizes serialization, criterion-view logic, outer-object trace schemas, lifecycle branch parity, and inner/joint trace branches. |

tests/families/test_general_family_mgcv_parity.py

@@ -53,6 +53,22 @@ def _gaulss_tensor_data(seed=22, n=160):
     return pd.DataFrame({"y": y, "x0": x0, "x1": x1})
 
 
+GENERAL_MULTISMOOTH_FORMULA = [
+    'y ~ s(x0, bs="cr", k=8) + s(x1, bs="cr", k=8)',
+    '~ s(x0, bs="cr", k=7) + s(x1, bs="cr", k=7)',
+]
+
+
+def _gaulss_multismooth_data(seed=1501, n=220):
+    rng = np.random.default_rng(seed)
+    x0 = rng.uniform(-1.5, 1.5, size=n)
+    x1 = rng.uniform(-1.5, 1.5, size=n)
+    mu = 0.3 + 1.8 * np.sin(2.5 * x0) + 1.1 * np.cos(3.0 * x1)
+    sigma = np.exp(-1.0 + 0.45 * np.sin(2.0 * x0) - 0.35 * np.cos(2.5 * x1))
+    y = rng.normal(mu, sigma, size=n)
+    return pd.DataFrame({"y": y, "x0": x0, "x1": x1})
+
+
 def _gammals_data(n=100, seed=2):
     rng = np.random.default_rng(seed)
     x = rng.uniform(-1.0, 1.0, n)

tests/optimization/test_mgcv_general_family_preoptimization_parity.py

@@ -20,10 +20,12 @@
 from tests._paths import PARITY_DIR, REPO_ROOT
 from tests.families.test_general_family_mgcv_parity import (
     GAULSS_FORMULA,
+    GENERAL_MULTISMOOTH_FORMULA,
     _gammals_by_data,
     _gammals_data,
     _gaulss_by_data,
     _gaulss_data,
+    _gaulss_multismooth_data,
 )
 from tests.mgcv_parity_utils import _family_specs, _fit_nampy_model_fixed_sp
 from tests.reference_fixtures import (
@@ -404,6 +406,15 @@ def _assert_general_fit5_setup_parity(
 
 GENERAL_PREOPT_CASES = [
     ("gaulss_cr", "gaulss", GAULSS_FORMULA, _gaulss_data, "ML", False, True),
+    (
+        "gaulss_multi_smooth_both_predictors",
+        "gaulss",
+        GENERAL_MULTISMOOTH_FORMULA,
+        _gaulss_multismooth_data,
+        "ML",
+        False,
+        True,
+    ),
     (
         "gaulss_fs",
         "gaulss",

tests/parity/test_mgcv_general_family_prediction_filters_parity.py

@@ -3,6 +3,7 @@
 import numpy as np
 import pytest
 
+from nampy.gam.inference.summary import summary_gam
 from tests.families.test_general_family_mgcv_parity import (
     _gaulss_by_data,
     _general_newdata,
@@ -17,6 +18,19 @@
 ]
 
 
+def test_general_family_inference_uses_predictor_aware_term_labels():
+    """Formula-list parametric and smooth rows use mgcv's later-LP suffix."""
+    data = _gaulss_by_data(seed=270, n=120)
+    gam = _fit_nampy_model(data, _FORMULA, "gaulss", "fixed")
+
+    summary = summary_gam(gam)
+    assert list(summary.pterms_table["label"]) == ["x", "z.1"]
+    assert list(summary.s_table["label"]) == ["s(z)", "s.1(x)"]
+    anova = gam.anova()
+    assert list(anova.parametric_table["label"]) == ["x", "z.1"]
+    assert list(anova.smooth_table["label"]) == ["s(z)", "s.1(x)"]
+
+
 def test_general_family_terms_filter_values_labels_and_se_match_mgcv():
     data = _gaulss_by_data(seed=271, n=120)
     newdata = _general_newdata(data, n=19)