Smooth terms

[BerriJ]

Implementation complete. All 871 tests pass (14 new), ruff clean. Summary of what was built:

Phase 1+3 — splines.py (new)

• make_bspline_knots / make_bspline_basis: equidistant-knot B-spline basis via scipy.interpolate.BSpline.design_matrix, with exact linear extrapolation beyond boundary knots (mgcv bs="ps" behavior)
• make_difference_penalty: Frobenius-normalized $S = D_q^\top D_q$
• make_centering_constraint: sum-to-zero identifiability via QR null space of the weighted basis column means ($K-1$ coefficients, no intercept clash)
• make_demmler_reinsch_transform: fixed-at-fit reparameterization giving identity Gramian + diagonal penalty
• Frozen PSplineTermState and PSplineTerm(Term) with feature, n_splines=20, degree=3, diff_order=2, lambda_=None (grid selection) or fixed, ic, forget, knot_padding=0.05; immutable update returning new instances; edf_ / lambda_selected_ diagnostics; EDF-based IC selection using $\mathrm{tr}((G+\lambda S)^{-1}G)$ and calculate_effective_training_length

Phase 2 — solver & method

• cd_base.py: numba online_coordinate_descent_quadratic (runs CD on $G + \lambda S$, keeps soft-thresholding for L1; reduces bit-for-bit to plain CD when $S=0$) and online_coordinate_descent_quadratic_path with get_start_beta warm starts
• quadratic_penalty.py: QuadraticPenaltyPath(EstimationMethod) — path-based, geometric λ-grid scaled by $\mathrm{tr}(G)/\mathrm{tr}(S)$ or fixed lambda_, reuses init_gram/update_gram with forgetting; registered as "quadratic_penalty" in the factory

Verification — test_terms_pspline.py: partition of unity, linear extrapolation, penalty annihilates degree-$(q{-}1)$ polynomials + bandedness, CD ≡ direct solve, $S=0$ ≡ plain CD (exact), warm-started path correctness, edf monotone in λ with bounds, large-λ polynomial limit, batch-vs-online Gram exactness, and full estimator integration on heteroskedastic $y \sim N(\sin 2x, e^{-0.5+0.3x})$ including update.

One notable deviation from the plan: plain CD on the raw centered basis is too ill-conditioned across the λ-path (descending warm starts hit max-iterations and return wrong solutions). I added the Demmler–Reinsch reparameterization fixed at fit time — Gramian becomes identity, penalty diagonal — so CD converges fast and exactly at every λ, while online Gram updates remain exact since the transform is a fixed linear map of the basis. Docs were added in terms_and_features.md and a runnable example in pspline_terms.py.