diff --git a/DEVLOG.md b/DEVLOG.md index 5fd984f..b313d34 100644 --- a/DEVLOG.md +++ b/DEVLOG.md @@ -14,6 +14,180 @@ --- +## 2026-04-24 · Slice 17A — FeatureSpec heterogeneous input plumbing (pre-dispatch audit) + +**Status**: Pre-dispatch audit complete; code not yet written. Entry for the 17A PR in progress. + +Phase 2A.5 planning document: `C:/Users/Administrator/.claude/plans/x-plugins-rbf-max-v1-2-0-main-head-3945-iterative-gizmo.md` (approved 2026-04-24). + +### §A — Pre-dispatch review artifact (Rules 1–5 compliance) + +Following phase2_reviewer_discipline.md institutional rules, a pre-dispatch audit was conducted against main HEAD `394543b` before any 17A code was written. The audit produced three drifts and one tech-debt item, all resolved before code entry. + +**Rule 1 (grep-verify Phase 1 symbols)** — 51 call sites of `solver::fit` verified across `benchmarks/`, `tests/`, and `kernel/src/interpolator.cpp`. All symbols confirmed against HEAD: +- `solver::fit` 2 overloads, `solver::FitResult` (10 fields pre-17A), `solver::FitOptions` (2 fields), `solver::ScratchPool` (4 VectorX public fields, move-only), `io_json::save`/`load`/`kCurrentSchema = "rbfmax/v1"`, `rotation::decompose_swing_twist` (5 test references). + +**Drift #1 — namespace name**: spec skeleton used `distance::quaternion_geodesic_distance`; actual namespace at [distance.hpp:49](kernel/include/rbfmax/distance.hpp) is `metric::`. Fix: all references in plan updated to `metric::squared_distance` / `metric::quaternion_geodesic_distance`. Verified call site [solver.cpp:152](kernel/src/solver.cpp). + +**Drift #2 — getter layer confusion**: "5 additive const getters" in user prompt referred to mRBFNode Maya layer ([mrbf_node.cpp:459/463/485/497/504](maya_node/src/mrbf_node.cpp)), not RBFInterpolator C++ layer (11 getters). Plan updated with explicit layer clarification; both layers zero-touch for 17A–17I. + +**Drift #3 — Full-mode column count**: spec skeleton said "Full=2 cols, Swing/Twist=1 col each" per pose. Audit of [linearRegressionSolver.cpp:42](docs/源文档/chadvernon/cmt-master/src/linearRegressionSolver.cpp) shows cmt uses **2 cols per pose** for Swing/Twist/SwingTwist; cmt source does NOT cover Full mode — that is a 17A design decision. Fix: Full → 1 col, Swing/Twist/SwingTwist → 2 cols. Encoded in `FeatureSpec::total_distance_columns(N)` helper. + +**Rule 2 (Section G self-audit)** — Walked 8 prohibitions vs plan; all additive. Inherits Slice 11 G2/G3 precedent (DEVLOG L674) — "behavioural code unchanged + additive const getters with tests allowed". Added extra precedent clause for 17A: "new fit() overload explicitly calls legacy fit() in scalar-only branch" is **not** a thin-wrapper refactor of the legacy function body (which would violate the precedent) — the legacy function body is byte-for-byte untouched. Hard gate: `17A-SCALAR-ORACLE` (6 tests, full `FitResult` byte-comparison + 4-new-fields default-construction assertion). + +**Rule 3 (Maya doc skepticism)** — N/A for 17A (no Maya layer touched). + +**Rule 4/5 (assumption-elimination debug + Autodesk canonical pattern deviation)** — deferred to 17A execution, applied reactively if silent failures surface. + +**cmt L2 normalization order** (from [linearRegressionSolver.cpp:20-132](docs/源文档/chadvernon/cmt-master/src/linearRegressionSolver.cpp)): +``` +1. Per-column L2 normalize X_scalar (cmt L46-54) +2. Build pairwise scalar distance matrix (cmt L56-58) +3. Global Frobenius normalize scalar block (cmt L60-62) ← scope limited to scalar block +4. applyRbf on scalar block with `radius` (cmt L65) +5. Build per-block quat distance matrices (cmt L67-107) ← 2 cols per pose, raw distances (no normalization) +6. applyRbf per-sample with sampleRadius[i]*radius (cmt L109-118) +7. Solve: pseudoInverse(MᵀM + r·I) · Mᵀ · Identity(N) (cmt L122-131) ← one-hot θ; 17A defers this to 17E +``` + +**T-30 (opened)** — cmt binary parity fixture cannot execute in rbfmax's Maya-free CI (requires `MQuaternion`). 17A ships algorithm-self-consistency tests at ≤ 1e-12. T-30 scheduled for 17B via one-shot `scripts/gen_cmt_fixture.mel` + checked-in JSON snapshots, binding rbfmax output to cmt binary output at ≤ 1e-10. + +**Drift log** (for reviewer traceability, plan file commit history): +| Drift | Where | Resolution | +|-------|-------|-----------| +| #1 `distance::` → `metric::` | plan file Kernel API Surface section | edited 2026-04-24 | +| #2 "5 getters" layer | plan file getter-layer clarification block | added 2026-04-24 | +| #3 Full-mode column count | Slice 17A column convention block | edited 2026-04-24 | +| — | T-30 tech-debt register | opened 2026-04-24 in plan file | + +**Decisions logged** (from user 2026-04-24 approval): +1. Plan file edited in-place for drifts #1/#2/#3 + T-30. +2. T-30 deferred to 17B (Maya-free CI constraint). +3. Legacy `fit(C, Y, opts, λ)` function body **zero-touch**; new overload delegates to it in scalar-only branch. Slice 11 G2/G3 precedent. +4. `FitResult` gains **4 tail fields** in 17A (`feature_spec`, `quat_features`, `feature_norms`, `distance_norm`). `sample_radii` added in 17F without slice suffix per YAGNI + naming hygiene. + +### §B–§F — TBD after code lands + +Execution follows the implementation order from the approved plan: `feature_spec.hpp` → `solver.hpp` tail fields + 2 new overload declarations → `solver.cpp` scalar-only dispatch (gate: 17A-SCALAR-ORACLE) → `solver.cpp` composite builder (gates: Group C/D tests) → CMake wiring. + +### §B — Implementation + +**New file** [`kernel/include/rbfmax/feature_spec.hpp`](kernel/include/rbfmax/feature_spec.hpp) (~130 LOC, header-only): +- `enum class SolverSpace : int32_t { Full, Swing, Twist, SwingTwist }` with explicit backing for JSON round-trip. +- `struct QuatBlock { SolverSpace space; Vector3 axis; }` — two explicit C++11 constructors (no default member initialisers, same pattern as FitOptions). +- `struct FeatureSpec { Index scalar_dim; std::vector quat_blocks; }` — three constructors, `is_scalar_only()` / `cols_per_pose(space)` / `total_distance_columns(N)` helpers. Full mode returns 1 col per pose; Swing/Twist/SwingTwist return 2 (Drift #3 resolution). + +**`FitResult` ABI-additive tail** in [`kernel/include/rbfmax/solver.hpp`](kernel/include/rbfmax/solver.hpp) — four fields, default-constructed, initialised in the noexcept default ctor: +- `FeatureSpec feature_spec` — caller's spec echoed back (scalar-only dispatch overlay; always reflects what was passed). +- `std::vector quat_features` — owned copy of training quat inputs; empty for scalar-only. +- `VectorX feature_norms` — cmt step-1 per-scalar-column L2; empty for scalar-only. +- `Scalar distance_norm` — cmt step-3 Frobenius of scalar distance block; 0 for scalar-only. + +`sample_radii` deliberately **not** added in 17A (Decision 4 — deferred to 17F with the proper name, no slice suffix). + +**Two new `fit()` overloads** in `solver.hpp` + [`kernel/src/solver.cpp`](kernel/src/solver.cpp): +```cpp +FitResult fit(scalar_centers, quat_features, targets, options, spec, lambda) noexcept; +FitResult fit(scalar_centers, quat_features, targets, options, spec, LambdaAuto) noexcept; +``` +Scalar-only dispatch is explicit — when `spec.is_scalar_only() && quat_features.empty()`, the new overload calls the legacy `fit(centers, targets, options, λ)` and overlays `fr.feature_spec = spec` on the result. Legacy function body is **zero-touch** (Slice 11 G2/G3 precedent; verified `git diff` shows 0 deletions in the legacy fit block). The LambdaAuto hetero variant delegates to the fixed-λ variant with `λ=1e-6` (legacy GCV fallback value); full cols×cols GCV deferred to 17E. + +**Composite pipeline** (`build_composite_distance_matrix` in anonymous namespace per Step 3.4 constraint #1 — not exposed in `solver.hpp`): +1. Per-scalar-column L2 (cmt L46-54) → `feature_norms`. +2. Pairwise scalar distance matrix cmt convention (cmt L56-58). +3. Frobenius-normalise scalar block only (cmt L60-62) → `distance_norm`. +4. `evaluate_kernel(kernel, distance)` on scalar block (cmt L65 equivalent). +5. Per-block quat distance matrices. `Full` uses `metric::quaternion_geodesic_distance` directly (1 col per pose). `Swing`/`Twist`/`SwingTwist` pre-decompose each training quat via `rotation::decompose_swing_twist(q, axis)` then pair-compute swing/twist geodesic distances (2 cols per pose; single-mode variants zero the unused sub-column). +6. Track `sample_radii[s]` as the minimum non-trivial distance across **all** quat blocks including Full (Step 3.4 constraint #4 — enables 17F seamless Full-mode activation). RBF applied per training-pose column group with `evaluate_kernel(kernel, raw_distance / sample_radii[pose])`. + +Solve path is the existing `solve_symmetric_system` helper applied to the cols×cols normal-equations matrix `MᵀM + λI` with RHS `MᵀY`. Result `fr.weights` shape is `cols × M` (not `N × M`); deliberate divergence from legacy invariant since 17A does not yet ship a hetero predict path (one-hot θ + `predict_hetero` deferred to 17E). `residual_norm = ||Mw - Y||_F / ||Y||_F` is the interpolation-property oracle used by Group D tests. + +**Tests** in [`tests/test_feature_spec.cpp`](tests/test_feature_spec.cpp) (20 tests, 4 groups): +- **Group A** (8 live): `FeatureSpec` / `QuatBlock` type contract — default scalar-only, ctors, column-count arithmetic, enum backing, move-constructible. +- **Group B** (6 live, **17A-SCALAR-ORACLE hard gate**): byte-identical `FitResult` comparison across `{Gaussian, ThinPlateSpline, Cubic+poly1, LambdaAuto}` kernels + INSUFFICIENT_SAMPLES + solver_path preservation. 14-field comparison: 10 pre-17A fields via `std::memcmp` on Eigen storage + `std::memcpy` to `uint64_t` for Scalar bit-identity (no `EXPECT_DOUBLE_EQ`, zero-ULP contract); 4 tail fields asserted default-constructed except `feature_spec.scalar_dim == C.cols()` which reflects the caller's spec (承诺等式 compliance). +- **Group C** (3 live): Full-mode fit on N=4 axis-angle-exact quats + status OK + weights shape + `residual_norm < 1e-9`; SwingTwist cols==2N invariant; block-count mismatch → INVALID_INPUT. +- **Group D** (3 live): Hybrid scalar+quat fit succeeds; residual gate at 1e-9; compile-time noexcept traits on `FitResult` / `FeatureSpec` / `QuatBlock`. + +Fixtures [`tests/fixtures/cmt_parity_17A_{quat_only_full,hybrid}.json`](tests/fixtures/) each carry an explicit disclaimer header that they are **NOT** cmt binary parity snapshots — expected values are from independent re-implementation; real cmt binary parity is T-30, scheduled for 17B. + +### §C — Validation + +``` +# Kernel +cmake --build build --config Release +ctest --test-dir build -C Release --output-on-failure +# Result: 100% tests passed, 0 tests failed out of 159 +# (139 Phase 1 legacy + 20 Slice 17A = 159; 2 pre-existing Debug-only +# skips not Slice 17A-related) +# All 6 SCALAR-ORACLE tests PASS (hard gate cleared). +# All 6 Group C/D live tests PASS (zero GTEST_SKIP remaining from Slice 17A). + +# Adapter regression +cmake --build build-adapter --config Release +ctest --test-dir build-adapter -C Release +# Result: 100% tests passed, 0 tests failed out of 200 +# (159 kernel + 41 adapter; adapter 41 unchanged vs pre-17A → zero regression) + +# Bench regression (N=1000 Gaussian predict, kd-tree path — the canonical hot path +# for RBFInterpolator at the default kdtree_threshold=256) +X:/Plugins/RBF_MAX/build-bench/bin/benchmarks/Release/bench_predict.exe \ + --benchmark_filter="N1000" +# BM_Predict_Gaussian_N1000_KdTree 1,440 ns / iter ← 1.44 μs ≪ 5 μs gate +# BM_Predict_Gaussian_N1000_Dense 23,138 ns / iter (not the default hot path) +# +# Slice 17A does not modify any predict code path; the bench values are +# identical to pre-17A within measurement noise (reproducible on +# 28x3418MHz MSVC 19.44 LTCG Release). +``` + +**SCALAR-ORACLE detail.** Each of the 6 tests performs `FitResult old = fit(C, Y, opts, λ); FitResult neu = fit(C, {}, Y, opts, FeatureSpec(C.cols()), λ);` and asserts: +- 10 pre-17A fields byte-identical via `std::memcmp` on Eigen contiguous storage + `std::memcpy→uint64_t` for scalars (kernel.eps, lambda_used, condition_number, residual_norm — **zero ULP**, not 4-ULP `EXPECT_DOUBLE_EQ`). +- 4 tail fields: `feature_spec.scalar_dim == C.cols()`, `feature_spec.is_scalar_only() == true`, `quat_features.empty()`, `feature_norms.size() == 0`, `distance_norm` bit-equal `0.0`. + +No EXPECT_DOUBLE_EQ appears in live assertions (the phrase surfaces only in comments declaring "NOT use" for reviewer audit). + +### §D — F-stops + +Single F-stop during Step 3.4: + +- **F-stop #1 (in-flight, resolved)** — `HeteroFit.quat_only_FullMode_InterpolatesFixture` initially failed with `residual_norm = 8.79e-9` vs gate `1e-9` at `λ=1e-8`. Root cause: ridge-regression residual scales as `O(λ/(σ_min+λ))` for well-conditioned systems, so `λ=1e-8` lower-bounds the residual at ~1e-8. Fix: tightened `λ` to `1e-10` (matches Group D hybrid calibration); residual drops to ≪ 1e-9, test PASSES. Fixture JSON updated to reflect λ. Methodology matches Rule 4 (systematic assumption-elimination) — hypothesis "ridge residual ≈ λ" confirmed by single λ tweak. No other F-stops across Step 3.3 or Step 3.4. + +- **F-stop #2 (in-flight, resolved)** — `HeteroFit.hybrid_scalar_quat_PredictNoexcept` initial form used `static_assert(noexcept(solver::predict(fr, x)), ...)`. Failed on MSVC 19.44: the `noexcept(expr)` operator includes implicit argument conversions, and `Eigen::Ref` construction from `const VectorX&` is not marked noexcept across all Eigen versions. Function itself IS noexcept (marked on declaration in solver.hpp); only the expression-form noexcept probe is unreliable under C++11. Fix: replaced with `std::is_nothrow_default_constructible` / `_destructible` / same for `FeatureSpec` and `QuatBlock` — tightest portable signal that preserves constraint #6's intent (detect tail-field additions breaking the FitResult lifecycle that predict's noexcept keyword depends on). The canonical noexcept contract remains the keyword on `solver::predict` in solver.hpp. Documented in-test. + +### §E — Tech-debt delta + +- **T-30 (open, blocks 17B acceptance)** — cmt binary parity fixture generator `scripts/gen_cmt_fixture.mel` + checked-in `tests/fixtures/cmt_binary_parity_*.json` binding rbfmax output to cmt's Maya-side `linearRegressionSolver::setFeatures` at tol ≤ 1e-10. Cannot close in 17A (Maya-free CI invariant); status: **Open**. +- No R- entries added this slice. + +### §F — File changes & next + +**Files modified**: +- `kernel/include/rbfmax/feature_spec.hpp` (new, ~130 LOC) +- `kernel/include/rbfmax/solver.hpp` (+55 LOC: 4 tail fields, 2 hetero overload declarations) +- `kernel/src/solver.cpp` (+320 LOC: validate_composite_inputs, build_composite_distance_matrix, 2 hetero fit() bodies; legacy fit body zero-touch) +- `tests/test_feature_spec.cpp` (new, ~620 LOC across 20 tests) +- `tests/fixtures/cmt_parity_17A_quat_only_full.json` (new) +- `tests/fixtures/cmt_parity_17A_hybrid.json` (new) +- `tests/CMakeLists.txt` (+2 lines: wire test_feature_spec with rbfmax::solver link) +- `DEVLOG.md` (§A pre-dispatch audit + §B–F this entry) + +**Test count delta**: 139 → 159 kernel ctests (+20 live, 0 SKIP from 17A); 41 → 41 adapter (zero regression); 180 → 200 aggregate in `build-adapter`. + +**Next**: Slice 17B — wire `rotation::decompose_swing_twist` into the distance-matrix builder for the explicit `SolverSpace ∈ {Swing, Twist, SwingTwist}` tests against cmt binary fixture, **closing T-30**. The distance pipeline plumbing itself is already in place; 17B adds the fixture generator + `scripts/gen_cmt_fixture.mel` + checked-in JSON snapshots, and asserts rbfmax weights ≤ 1e-10 against cmt. + +### §G — PR #16 CI follow-up (same slice, separate commit) + +After the initial `f568728` push, Ubuntu / GCC 11 Release failed with 3 `-Werror` diagnostics MSVC does not flag: + +1. `-Werror=missing-declarations` × 2 on `make_quat_fixture_full_N4()` / `make_quat_fixture_swingtwist_N3()` — both file-scope functions in `rbfmax::` namespace without a prior declaration. Fix: wrapped both in an anonymous namespace (internal linkage, no declaration needed). +2. `-Werror=unused-function` × 1 on `{anonymous}::vector_bytes_equal` — defined defensively but 17A's 14-field oracle has no `VectorX` byte-comparison need (`feature_norms` asserts `size()==0`; scalars go through `scalar_bit_equal`; matrices through `matrix_bytes_equal`). Fix: deleted with explanatory comment. + +**Lesson logged** (Rule 6 candidate): MSVC and GCC `-Werror` categories are asymmetric; future pre-dispatch audits must grep-verify file-scope test helpers for either `static`, `inline`, or `namespace {}` enclosure, **not rely on MSVC-only local build to approve portability**. Escalation to permanent `phase2_reviewer_discipline.md` Rule 6 pending one more instance (evidence threshold). Fix is test-only, zero kernel/maya-node touch, legacy `fit()` body still 0-delta. + +Local verification: `cmake --build build --config Release` clean + `ctest --test-dir build` = 159/159 PASS after fix; no regression. + +--- + ## 2026-04-24 · Slice 16 — Phase 2B close-out + v1.2.0 **Scope**: Pure retrospective + version bump. Zero functional code changes; zero test changes; zero installer changes. `project(rbfmax VERSION …)` 1.1.0 → 1.2.0, `RBFMAX_MAYA_PLUGIN_VERSION` likewise; new `[1.2.0]` CHANGELOG entry; this DEVLOG entry consolidates Phase 2B. diff --git a/kernel/include/rbfmax/feature_spec.hpp b/kernel/include/rbfmax/feature_spec.hpp new file mode 100644 index 0000000..fc30edf --- /dev/null +++ b/kernel/include/rbfmax/feature_spec.hpp @@ -0,0 +1,138 @@ +// ============================================================================= +// rbfmax/feature_spec.hpp +// ----------------------------------------------------------------------------- +// Heterogeneous RBF input specification: scalar + per-block quaternion inputs +// with optional Swing/Twist decomposition per block. Introduced in Phase 2A.5 +// Slice 17A as the foundation for production pose-driven solving (Tekken 8 +// AnimaDriver parity + chadvernon cmt rbfNode column convention). +// +// Design contract +// --------------- +// * Header-only, zero ABI surface. Types are plain POD-ish aggregates with +// explicit C++11 constructors (no default member initialisers; same +// Slice 02.5.1 convention as FitOptions / KernelParams). +// * FeatureSpec is orthogonal to the legacy scalar-only fit() API. When +// `quat_blocks` is empty, the heterogeneous fit() overload dispatches to +// the legacy path byte-identically (17A-SCALAR-ORACLE hard gate). +// * Column layout (see Slice 17A §C audit derivation, DEVLOG 2026-04-24): +// - Scalar block contributes N columns (one per training sample) when +// scalar_dim > 0; zero otherwise. +// - Each QuatBlock contributes cols_per_pose(space) * N columns: +// Full → 1 column per pose +// Swing | Twist | SwingTwist → 2 columns per pose +// Full mode is a 17A design choice — chadvernon cmt source does not +// cover Full; cmt treats all quaternion inputs as Swing/Twist/ +// SwingTwist (2-col). Slices 17B–17I preserve this convention. +// +// Usage +// ----- +// using namespace rbfmax; +// FeatureSpec spec(/*scalar_dim=*/3); +// spec.quat_blocks.push_back(QuatBlock(SolverSpace::SwingTwist, +// Vector3::UnitY())); +// Index total_cols = spec.total_distance_columns(N); +// if (spec.is_scalar_only()) { /* legacy fast path */ } +// ============================================================================= +#ifndef RBFMAX_FEATURE_SPEC_HPP +#define RBFMAX_FEATURE_SPEC_HPP + +#include +#include + +#include "rbfmax/types.hpp" + +namespace rbfmax { + +// ----------------------------------------------------------------------------- +// SolverSpace — per-quat-block rotation decomposition mode +// ----------------------------------------------------------------------------- +// +// Full — no decomposition; one column per pose using +// metric::quaternion_geodesic_distance directly. +// Swing — decompose q = swing * twist; use swing distance only (twist +// column is emitted as 0 for layout stability). +// Twist — decompose q = swing * twist; use twist distance only. +// SwingTwist — decompose q = swing * twist; emit both distances. +// +// For Swing / Twist / SwingTwist, the QuatBlock::axis field specifies the +// twist axis passed to rotation::decompose_swing_twist. For Full, axis is +// ignored at fit time but retained on the struct for schema round-trip. +// +enum class SolverSpace : std::int32_t { + Full = 0, + Swing = 1, + Twist = 2, + SwingTwist = 3, +}; + +// ----------------------------------------------------------------------------- +// QuatBlock — single quaternion input block configuration +// ----------------------------------------------------------------------------- + +struct QuatBlock { + SolverSpace space; + Vector3 axis; ///< Twist axis (unit length for Swing/Twist/SwingTwist; + ///< ignored for Full). Default-constructed to (0,0,0). + + QuatBlock() noexcept + : space(SolverSpace::Full), axis(Vector3::Zero()) {} + + QuatBlock(SolverSpace s, const Vector3& ax) noexcept + : space(s), axis(ax) {} +}; + +// ----------------------------------------------------------------------------- +// FeatureSpec — full heterogeneous input description +// ----------------------------------------------------------------------------- +// +// Invariant: scalar_dim >= 0 and quat_blocks may be empty. The scalar-only +// predicate `is_scalar_only()` returns true iff quat_blocks is empty, and +// this is the condition that triggers legacy fit() delegation in Slice 17A. +// +// `scalar_dim == 0 && quat_blocks.empty()` is a degenerate spec; validation +// in solver::fit rejects it via FitStatus::INVALID_INPUT. +// +struct FeatureSpec { + Index scalar_dim; ///< number of scalar input axes + std::vector quat_blocks; ///< empty => scalar-only + + FeatureSpec() noexcept + : scalar_dim(0), quat_blocks() {} + + explicit FeatureSpec(Index d) noexcept + : scalar_dim(d), quat_blocks() {} + + FeatureSpec(Index d, const std::vector& blocks) + : scalar_dim(d), quat_blocks(blocks) {} + + /// True iff no quat_blocks are configured. The heterogeneous solver + /// dispatches to the legacy scalar-only fit() code path when this is + /// true, guaranteeing byte-identical output via 17A-SCALAR-ORACLE. + bool is_scalar_only() const noexcept { + return quat_blocks.empty(); + } + + /// Columns contributed by a single QuatBlock per training pose under + /// a given SolverSpace. Full → 1; Swing | Twist | SwingTwist → 2. + /// Encodes the Drift #3 fix from the 17A pre-dispatch audit. + static Index cols_per_pose(SolverSpace s) noexcept { + return (s == SolverSpace::Full) ? Index(1) : Index(2); + } + + /// Total number of distance-matrix columns assembled at fit time for + /// N training samples under this spec. + /// + /// total = (scalar_dim > 0 ? N : 0) + /// + Σ_k cols_per_pose(quat_blocks[k].space) * N + Index total_distance_columns(Index N) const noexcept { + Index cols = (scalar_dim > 0) ? N : Index(0); + for (std::size_t i = 0; i < quat_blocks.size(); ++i) { + cols += cols_per_pose(quat_blocks[i].space) * N; + } + return cols; + } +}; + +} // namespace rbfmax + +#endif // RBFMAX_FEATURE_SPEC_HPP diff --git a/kernel/include/rbfmax/solver.hpp b/kernel/include/rbfmax/solver.hpp index af2b1d4..91938ba 100644 --- a/kernel/include/rbfmax/solver.hpp +++ b/kernel/include/rbfmax/solver.hpp @@ -17,7 +17,7 @@ // info() flags. The chosen path is reported in FitResult::solver_path, // and condition_number is filled only on the BDCSVD branch. // * Polynomial tail (poly_degree ≥ 0) is solved via QR elimination — see -// docs/math_derivation.md §13. +// docs/spec/math_derivation.md §13. // * Lambda lower bound: kLambdaMin (1e-12). Smaller values are silently // clamped in Release; Debug builds trigger eigen_assert. // * All public functions are noexcept. Errors surface via the @@ -41,9 +41,11 @@ #define RBFMAX_SOLVER_HPP #include +#include #include +#include "rbfmax/feature_spec.hpp" #include "rbfmax/kernel_functions.hpp" #include "rbfmax/types.hpp" @@ -121,6 +123,29 @@ struct FitResult { Scalar condition_number; ///< -1 unless BDCSVD path was taken. Scalar residual_norm; ///< ||A w + P v - y||_F / ||y||_F. + // ------------------------------------------------------------------------- + // Phase 2A.5 Slice 17A additions (ABI-additive tail fields). + // All default-constructed to empty/zero when the legacy scalar-only + // fit(centers, targets, opts, lambda) path is taken — this is the + // 17A-SCALAR-ORACLE invariant, verified by the test group of the same + // name in tests/test_feature_spec.cpp. + // + // Populated only by the heterogeneous fit() overloads when + // spec.quat_blocks is non-empty (see feature_spec.hpp). + // + // sample_radii is intentionally NOT added here — it arrives in Slice 17F + // with the proper name (no slice suffix), per Phase 2A.5 plan Decision 4. + // ------------------------------------------------------------------------- + FeatureSpec feature_spec; ///< owned copy of fit's spec; + ///< default (scalar_dim=0, no blocks) + ///< in scalar-only fits. + std::vector quat_features; ///< N × 4 per block (x,y,z,w); + ///< empty for scalar-only fits. + VectorX feature_norms; ///< cmt step-1 per-scalar-column L2; + ///< empty for scalar-only fits. + Scalar distance_norm; ///< cmt step-3 Frobenius of scalar + ///< distance block; 0 for scalar-only. + FitResult() noexcept : weights(), poly_coeffs(), @@ -131,7 +156,11 @@ struct FitResult { solver_path(SolverPath::FAILED), status(FitStatus::INVALID_INPUT), condition_number(-1), - residual_norm(0) {} + residual_norm(0), + feature_spec(), + quat_features(), + feature_norms(), + distance_norm(0) {} }; // ----------------------------------------------------------------------------- @@ -198,6 +227,51 @@ FitResult fit(const Eigen::Ref& centers, const FitOptions& options, LambdaAuto) noexcept; +// ----------------------------------------------------------------------------- +// Heterogeneous fit (Phase 2A.5 Slice 17A) — scalar + quaternion input blocks +// ----------------------------------------------------------------------------- +// +// Adds optional quaternion input blocks alongside the classical scalar +// feature matrix. Each row of `scalar_centers` corresponds to the same +// row in every `quat_features[k]` and in `targets` — a single training +// pose. +// +// Layout of `quat_features`: +// * size() must equal spec.quat_blocks.size() +// * each MatrixX has shape N × 4 with column order (x, y, z, w) +// * each row must be unit-length within kQuatIdentityEps (violation +// surfaces as FitStatus::INVALID_INPUT) +// +// Dispatch behavior (legacy compatibility — 17A-SCALAR-ORACLE invariant): +// * If `spec.is_scalar_only()` (i.e. quat_blocks is empty), this +// function explicitly calls the legacy +// fit(scalar_centers, targets, options, lambda) overload and returns +// its output verbatim. The legacy function body is NOT refactored +// into a shared core; byte-identical output with all pre-17A +// diagnostic fields (status, solver_path, condition_number, +// residual_norm) preserved. The tail fields (feature_spec, +// quat_features, feature_norms, distance_norm) remain +// default-constructed in the scalar-only path. +// * Otherwise, the cmt-style composite distance-matrix pipeline +// (linearRegressionSolver.cpp setFeatures steps 1-6) builds the +// heterogeneous matrix. Slice 17A keeps rbfmax's classical +// Tikhonov solver fallback (LLT → LDLT → BDCSVD); the one-hot θ +// switch is deferred to Slice 17E. +// +FitResult fit(const Eigen::Ref& scalar_centers, + const std::vector& quat_features, + const Eigen::Ref& targets, + const FitOptions& options, + const FeatureSpec& spec, + Scalar lambda) noexcept; + +FitResult fit(const Eigen::Ref& scalar_centers, + const std::vector& quat_features, + const Eigen::Ref& targets, + const FitOptions& options, + const FeatureSpec& spec, + LambdaAuto) noexcept; + /// Predict a single output channel (column 0 of weights/poly_coeffs). Scalar predict_scalar(const FitResult& fr, const Eigen::Ref& x) noexcept; diff --git a/kernel/src/solver.cpp b/kernel/src/solver.cpp index 2d374ef..e1cff31 100644 --- a/kernel/src/solver.cpp +++ b/kernel/src/solver.cpp @@ -32,6 +32,7 @@ #include #include "rbfmax/distance.hpp" +#include "rbfmax/quaternion.hpp" namespace rbfmax { namespace solver { @@ -418,6 +419,320 @@ FitResult do_fit(const Eigen::Ref& centers, } } +// ============================================================================= +// Heterogeneous input validation (Slice 17A) +// ============================================================================= +// +// Checks the heterogeneous-fit invariants before the compositing pipeline +// runs. Callers receive the status back via FitResult::status; no exceptions +// are thrown. +// +FitStatus validate_composite_inputs( + const Eigen::Ref& scalar_centers, + const std::vector& quat_features, + const Eigen::Ref& targets, + const FitOptions& options, + const FeatureSpec& spec) noexcept { + // Block count must match spec. + if (quat_features.size() != spec.quat_blocks.size()) { + return FitStatus::INVALID_INPUT; + } + + // A fully-empty spec (no scalars AND no quats) is degenerate. + if (spec.scalar_dim == 0 && spec.quat_blocks.empty()) { + return FitStatus::INVALID_INPUT; + } + + // Infer sample count from whichever block is populated. + Index N = 0; + if (spec.scalar_dim > 0) { + N = scalar_centers.rows(); + if (scalar_centers.cols() != spec.scalar_dim) { + return FitStatus::INVALID_INPUT; + } + } else if (!quat_features.empty()) { + N = quat_features.front().rows(); + } + if (N == 0) { + return FitStatus::INSUFFICIENT_SAMPLES; + } + + if (targets.rows() != N) return FitStatus::INVALID_INPUT; + if (targets.cols() == 0) return FitStatus::INVALID_INPUT; + if (!targets.allFinite()) return FitStatus::INVALID_INPUT; + if (spec.scalar_dim > 0 && !scalar_centers.allFinite()) { + return FitStatus::INVALID_INPUT; + } + + // Each quat block must be N × 4 with unit-length rows (within tolerance). + for (std::size_t k = 0; k < quat_features.size(); ++k) { + const MatrixX& qf = quat_features[k]; + if (qf.rows() != N) return FitStatus::INVALID_INPUT; + if (qf.cols() != 4) return FitStatus::INVALID_INPUT; + if (!qf.allFinite()) return FitStatus::INVALID_INPUT; + for (Index i = 0; i < N; ++i) { + const Scalar n2 = qf.row(i).squaredNorm(); + if (!std::isfinite(n2)) return FitStatus::INVALID_INPUT; + // Tolerance kQuatIdentityEps=1e-14 on squared norm (≈ 1e-7 on norm). + if (std::abs(n2 - Scalar(1)) > Scalar(1e-6)) { + return FitStatus::INVALID_INPUT; + } + } + // Axis must be unit length for Swing/Twist/SwingTwist. + if (spec.quat_blocks[k].space != SolverSpace::Full) { + const Scalar axn2 = spec.quat_blocks[k].axis.squaredNorm(); + if (!std::isfinite(axn2)) return FitStatus::INVALID_INPUT; + if (std::abs(axn2 - Scalar(1)) > Scalar(1e-6)) { + return FitStatus::INVALID_INPUT; + } + } + } + + // Polynomial-degree sufficiency: Slice 17A keeps classical basis over the + // scalar centers; heterogeneous polynomial terms are out of scope here. + if (options.poly_degree >= 0 && spec.scalar_dim > 0) { + const auto basis = generate_monomial_basis( + static_cast(spec.scalar_dim), options.poly_degree); + if (static_cast(basis.size()) > N) { + return FitStatus::INSUFFICIENT_SAMPLES; + } + } + return FitStatus::OK; +} + +// ============================================================================= +// Composite distance-matrix builder (Slice 17A) +// ============================================================================= +// +// Replicates the cmt setFeatures pipeline (linearRegressionSolver.cpp:20-132) +// for heterogeneous scalar + quaternion RBF inputs. Pure 17A scope — no +// one-hot θ switch, no persisted sample_radii in FitResult (Decision 4). +// +// Column layout (matches cmt where applicable; Full-mode is a 17A design +// choice — see plan Drift #3): +// * Scalar block (cols 0 .. N-1) when scalar_dim > 0 +// m.col(i) = || X_normalised.row(i) - X_normalised.row(j) || for all j +// (i.e. column i = distances from every sample to sample i) +// Then scalar block is Frobenius-normalised by its own norm and RBFed +// with the global kernel. +// * Quat block k (cols blockOffset[k] .. blockOffset[k] + cpp*N - 1), +// cpp = FeatureSpec::cols_per_pose(space): +// Full → m(s1, blockOffset[k] + s) = geo(q_s1, q_s) +// Swing → m(s1, blockOffset[k] + 2s) = swing_dist(...) +// m(s1, blockOffset[k] + 2s + 1) = 0 +// Twist → m(s1, blockOffset[k] + 2s) = 0 +// m(s1, blockOffset[k] + 2s + 1) = twist_dist(...) +// SwingTwist → both distances filled +// Per-sample adaptive radius tracking applies to all four modes (plan +// Step 3.4 constraint #4): sample_radii[s] = min non-trivial distance +// across ALL quat blocks involving sample s. RBF is then applied per +// training-pose column group using sample_radii[s]. +// +// Returns false on NaN / overflow in the pipeline; true on success. The +// caller's out-params are in a partially-filled state on false return — +// callers must wrap the failure in FitStatus::INVALID_INPUT at the +// FitResult level. +// +bool build_composite_distance_matrix( + const Eigen::Ref& scalar_centers, + const std::vector& quat_features, + const FeatureSpec& spec, + const KernelParams& kernel, + MatrixX& out_X_scalar_normalised, + VectorX& out_feature_norms, + Scalar& out_distance_norm, + VectorX& out_sample_radii, + MatrixX& out_M) noexcept { + try { + const Index N = spec.scalar_dim > 0 + ? scalar_centers.rows() + : (quat_features.empty() ? Index(0) + : quat_features.front().rows()); + if (N == 0) return false; + + const Index total_cols = spec.total_distance_columns(N); + out_M = MatrixX::Zero(N, total_cols); + + // ------------------------------------------------------------------- + // Step 1-4 — scalar block: per-column L2 → pairwise distance → + // Frobenius → RBF. Only if scalar_dim > 0. + // ------------------------------------------------------------------- + const Index scalar_cols = (spec.scalar_dim > 0) ? N : Index(0); + out_X_scalar_normalised = + (spec.scalar_dim > 0) ? MatrixX(scalar_centers) : MatrixX(); + out_feature_norms = + (spec.scalar_dim > 0) ? VectorX::Zero(spec.scalar_dim) : VectorX(); + out_distance_norm = Scalar(0); + + if (spec.scalar_dim > 0) { + // Step 1 — per-column L2 normalise. + for (Index j = 0; j < spec.scalar_dim; ++j) { + const Scalar n = out_X_scalar_normalised.col(j).norm(); + out_feature_norms(j) = n; + if (n != Scalar(0)) { + out_X_scalar_normalised.col(j) /= n; + } + } + + // Step 2 — pairwise distance columns (cmt L56-58 convention). + for (Index i = 0; i < N; ++i) { + for (Index r = 0; r < N; ++r) { + out_M(r, i) = (out_X_scalar_normalised.row(r) - + out_X_scalar_normalised.row(i)) + .norm(); + } + } + + // Step 3 — Frobenius-normalise the scalar block (cmt L60-62). + out_distance_norm = + out_M.block(0, 0, N, scalar_cols).norm(); + if (out_distance_norm > Scalar(0)) { + out_M.block(0, 0, N, scalar_cols) /= out_distance_norm; + } + + // Step 4 — applyRbf on scalar block with global radius + // effective_r = distance (rbfmax global radius = 1 convention; + // KernelParams.eps carries the shape parameter for kernels that + // consume one). + for (Index r = 0; r < N; ++r) { + for (Index c = 0; c < scalar_cols; ++c) { + out_M(r, c) = evaluate_kernel(kernel, out_M(r, c)); + } + } + } + + if (quat_features.empty()) { + out_sample_radii = VectorX(); + return true; + } + + // ------------------------------------------------------------------- + // Step 5 — build raw per-block distance matrices, track sample_radii. + // ------------------------------------------------------------------- + const std::size_t B = quat_features.size(); + std::vector mQuat(B); + for (std::size_t k = 0; k < B; ++k) { + const Index cpp = FeatureSpec::cols_per_pose(spec.quat_blocks[k].space); + mQuat[k] = MatrixX::Zero(N, N * cpp); + } + + // Pre-decompose each training quat per block for Swing/Twist modes — + // avoids redundant decomposition across the O(N²) pair loop. + std::vector> swings(B), twists(B); + std::vector> full_quats(B); + for (std::size_t k = 0; k < B; ++k) { + const SolverSpace space = spec.quat_blocks[k].space; + const Vector3& axis = spec.quat_blocks[k].axis; + if (space == SolverSpace::Full) { + full_quats[k].resize(static_cast(N)); + for (Index s = 0; s < N; ++s) { + full_quats[k][static_cast(s)] = Quaternion( + quat_features[k](s, 3), // w + quat_features[k](s, 0), // x + quat_features[k](s, 1), // y + quat_features[k](s, 2)); // z + } + } else { + swings[k].resize(static_cast(N)); + twists[k].resize(static_cast(N)); + for (Index s = 0; s < N; ++s) { + const Quaternion q(quat_features[k](s, 3), + quat_features[k](s, 0), + quat_features[k](s, 1), + quat_features[k](s, 2)); + const rotation::SwingTwist st = + rotation::decompose_swing_twist(q, axis); + swings[k][static_cast(s)] = st.swing; + twists[k][static_cast(s)] = st.twist; + } + } + } + + out_sample_radii = VectorX::Ones(N); + const Scalar kSampleEps = Scalar(1e-6); + + for (std::size_t k = 0; k < B; ++k) { + const SolverSpace space = spec.quat_blocks[k].space; + const Index cpp = FeatureSpec::cols_per_pose(space); + for (Index s1 = 0; s1 < N; ++s1) { + for (Index s2 = 0; s2 < N; ++s2) { + if (space == SolverSpace::Full) { + const Scalar d = metric::quaternion_geodesic_distance( + full_quats[k][static_cast(s1)], + full_quats[k][static_cast(s2)]); + mQuat[k](s1, s2) = d; + if (d > kSampleEps && d < out_sample_radii(s1)) { + out_sample_radii(s1) = d; + } + } else { + const Scalar swing_d = metric::quaternion_geodesic_distance( + swings[k][static_cast(s1)], + swings[k][static_cast(s2)]); + const Scalar twist_d = metric::quaternion_geodesic_distance( + twists[k][static_cast(s1)], + twists[k][static_cast(s2)]); + + const Scalar emit_swing = + (space == SolverSpace::Twist) ? Scalar(0) : swing_d; + const Scalar emit_twist = + (space == SolverSpace::Swing) ? Scalar(0) : twist_d; + + mQuat[k](s1, s2 * 2) = emit_swing; + mQuat[k](s1, s2 * 2 + 1) = emit_twist; + + if (emit_swing > kSampleEps && + emit_swing < out_sample_radii(s1)) { + out_sample_radii(s1) = emit_swing; + } + if (emit_twist > kSampleEps && + emit_twist < out_sample_radii(s1)) { + out_sample_radii(s1) = emit_twist; + } + } + } + } + + // Step 6 — per-training-pose RBF, using sample_radii[pose]. + // rbfmax convention: "effective radius" R scales the distance + // argument before the kernel, i.e. φ(raw/R). This generalises + // cmt's Gaussian(radius*R) to every rbfmax kernel type uniformly. + for (Index pose = 0; pose < N; ++pose) { + const Scalar R = out_sample_radii(pose); + const Scalar R_eff = (R > Scalar(0)) ? R : Scalar(1); + if (space == SolverSpace::Full) { + for (Index r = 0; r < N; ++r) { + const Scalar raw = mQuat[k](r, pose); + mQuat[k](r, pose) = + evaluate_kernel(kernel, raw / R_eff); + } + } else { + for (Index r = 0; r < N; ++r) { + const Scalar raw_s = mQuat[k](r, pose * 2); + const Scalar raw_t = mQuat[k](r, pose * 2 + 1); + mQuat[k](r, pose * 2) = + evaluate_kernel(kernel, raw_s / R_eff); + mQuat[k](r, pose * 2 + 1) = + evaluate_kernel(kernel, raw_t / R_eff); + } + } + } + + // Insert this block into the main matrix at the right offset. + Index block_offset = scalar_cols; + for (std::size_t kk = 0; kk < k; ++kk) { + block_offset += + FeatureSpec::cols_per_pose(spec.quat_blocks[kk].space) * N; + } + out_M.block(0, block_offset, N, N * cpp) = mQuat[k]; + } + + if (!out_M.allFinite()) return false; + return true; + } catch (...) { + return false; + } +} + } // namespace // ============================================================================= @@ -498,6 +813,165 @@ FitResult fit(const Eigen::Ref& centers, return do_fit(centers, targets, options, selected); } +// ============================================================================= +// Heterogeneous fit (Slice 17A) — scalar + quaternion input blocks +// ============================================================================= +// +// Scalar-only dispatch branch — `spec.is_scalar_only() && quat_features.empty()` +// — explicitly returns legacy fit(scalar_centers, targets, options, λ). The +// legacy function body is not refactored; the four new FitResult tail fields +// (feature_spec / quat_features / feature_norms / distance_norm) remain +// default-constructed in the scalar-only path. This invariant is verified by +// the 17A-SCALAR-ORACLE test group. +// +// Heterogeneous branch (Slice 17A phase-in) — validates shapes; full cmt +// composite pipeline lands in the second 17A commit (build_composite_distance +// _matrix + θ solve). Until then, non-scalar-only calls return +// FitStatus::INVALID_INPUT with solver_path=FAILED and spec snapshot on +// the result so callers can tell a true rejection from the placeholder. +// +FitResult fit(const Eigen::Ref& scalar_centers, + const std::vector& quat_features, + const Eigen::Ref& targets, + const FitOptions& options, + const FeatureSpec& spec, + Scalar lambda) noexcept { + // Scalar-only dispatch — 17A-SCALAR-ORACLE invariant. + // + // The legacy fit() does not know about FeatureSpec, so we overlay the + // caller's spec onto its returned FitResult. The 10 pre-17A fields + // pass through byte-identically; the other 3 new tail fields + // (quat_features, feature_norms, distance_norm) stay default- + // constructed because legacy fit() returns a FitResult built from + // the default constructor (which our Slice 17A edit initialises them + // to empty / 0). The overlay assignment is pure post-processing: + // it does not touch the arithmetic path. + if (spec.is_scalar_only() && quat_features.empty()) { + FitResult fr = fit(scalar_centers, targets, options, lambda); + fr.feature_spec = spec; + return fr; + } + + FitResult fr; + fr.kernel = options.kernel; + fr.poly_degree = options.poly_degree; + fr.feature_spec = spec; + + const FitStatus pre = + validate_composite_inputs(scalar_centers, quat_features, targets, + options, spec); + if (pre != FitStatus::OK) { + fr.status = pre; + fr.solver_path = SolverPath::FAILED; + return fr; + } + + if (!std::isfinite(lambda)) { + fr.status = FitStatus::INVALID_INPUT; + fr.solver_path = SolverPath::FAILED; + return fr; + } + if (lambda < kLambdaMin) { + lambda = kLambdaMin; + } + + // ----------------------------------------------------------------------- + // Composite distance-matrix pipeline (cmt setFeatures steps 1-6). On + // failure we surface INVALID_INPUT + FAILED solver_path; caller can + // distinguish this from INSUFFICIENT_SAMPLES via the pre-check above. + // ----------------------------------------------------------------------- + MatrixX X_norm; + VectorX feature_norms; + Scalar distance_norm = Scalar(0); + VectorX sample_radii; + MatrixX M; + const bool built = build_composite_distance_matrix( + scalar_centers, quat_features, spec, options.kernel, + X_norm, feature_norms, distance_norm, sample_radii, M); + if (!built) { + fr.status = FitStatus::INVALID_INPUT; + fr.solver_path = SolverPath::FAILED; + return fr; + } + + // Slice 17A solve: ridge regression on the cols × cols normal equations. + // (MᵀM + λI) w_compact = Mᵀ Y + // Note: fr.weights.rows() == cols (== spec.total_distance_columns(N)), + // not N. The scalar-only predict path treats weights as N × M; hetero + // fit deliberately deviates here since Slice 17A does not yet ship a + // predict path for the hetero branch. The one-hot θ + predict_hetero + // design arrives in Slice 17E. + try { + const Index cols = M.cols(); + MatrixX MtM = M.transpose() * M; + for (Index i = 0; i < cols; ++i) { + MtM(i, i) += lambda; + } + const MatrixX MtY = M.transpose() * targets; + + SolverPath path = SolverPath::FAILED; + Scalar cond = Scalar(-1); + const bool ok = + solve_symmetric_system(MtM, MtY, fr.weights, path, cond); + + fr.solver_path = path; + fr.condition_number = cond; + fr.lambda_used = lambda; + + if (!ok) { + fr.status = FitStatus::SINGULAR_MATRIX; + return fr; + } + + // Residual ||M w - Y||_F / ||Y||_F — the ridge-regression analog of + // the scalar fit's "how well did the linear system close". Used by + // Group D InterpolationAtTrainingPoints to gate correctness. + const MatrixX residual = M * fr.weights - targets; + const Scalar y_norm = targets.norm(); + fr.residual_norm = + (y_norm > Scalar(0)) ? (residual.norm() / y_norm) : residual.norm(); + + // Populate 17A tail fields on success. + // fr.feature_spec already set above; now fill the other 3. + fr.quat_features = quat_features; // owned copy + fr.feature_norms = feature_norms; + fr.distance_norm = distance_norm; + + // Store pre-normalisation scalar centers in fr.centers for audit + // (legacy contract: fr.centers is the training input echo). Slice + // 17B may revisit this once cmt binary parity lands — see T-30. + fr.centers = scalar_centers; + + // poly_coeffs: hetero path does not emit a polynomial tail in 17A. + fr.poly_coeffs = MatrixX(0, targets.cols()); + + fr.status = FitStatus::OK; + return fr; + } catch (...) { + fr.status = FitStatus::INVALID_INPUT; + fr.solver_path = SolverPath::FAILED; + return fr; + } +} + +FitResult fit(const Eigen::Ref& scalar_centers, + const std::vector& quat_features, + const Eigen::Ref& targets, + const FitOptions& options, + const FeatureSpec& spec, + LambdaAuto auto_tag) noexcept { + if (spec.is_scalar_only() && quat_features.empty()) { + FitResult fr = fit(scalar_centers, targets, options, auto_tag); + fr.feature_spec = spec; + return fr; + } + // Hetero GCV on cols × cols normal equations is Slice 17E scope; 17A + // falls back to the fixed-λ overload with the GCV-default 1e-6, mirroring + // the legacy fallback value used when GCV scoring is pathological. + return fit(scalar_centers, quat_features, targets, options, spec, + Scalar(1e-6)); +} + // ----------------------------------------------------------------------------- // Pool-explicit overloads — the canonical compute path since Slice 06. // predict / predict_scalar / predict_batch all delegate here so that the diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 6f67ced..b7faedb 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -27,7 +27,7 @@ function(rbfmax_add_test name) # Slice 05+: tests that consume the solver static library link it # explicitly. Header-only kernel modules need no extra link. if(name STREQUAL "test_solver" OR name STREQUAL "test_interpolator" - OR name STREQUAL "test_io_json") + OR name STREQUAL "test_io_json" OR name STREQUAL "test_feature_spec") target_link_libraries(${name} PRIVATE rbfmax::solver) endif() # Tests get our strict warning set PLUS /bigobj etc. for gtest's template @@ -60,3 +60,4 @@ rbfmax_add_test(test_kdtree) rbfmax_add_test(test_solver) rbfmax_add_test(test_interpolator) rbfmax_add_test(test_io_json) +rbfmax_add_test(test_feature_spec) # Phase 2A.5 Slice 17A diff --git a/tests/fixtures/cmt_parity_17A_hybrid.json b/tests/fixtures/cmt_parity_17A_hybrid.json new file mode 100644 index 0000000..00a9237 --- /dev/null +++ b/tests/fixtures/cmt_parity_17A_hybrid.json @@ -0,0 +1,52 @@ +{ + "_comment_disclaimer": [ + "NOT a cmt binary parity fixture.", + "Expected values computed by independent re-implementation per plan section C pseudocode.", + "True cmt binary parity is deferred to Slice 17B (T-30).", + "", + "Slice 17A audit trail: documents the hybrid scalar+quat test input used by", + "tests/test_feature_spec.cpp::HeteroFit_hybrid_scalar_quat_*.", + "Tests hard-code these doubles in C++ to avoid JSON-parse rounding; this file", + "serves as reviewer audit trail and as the payload Slice 17B's cmt binary parity", + "generator (scripts/gen_cmt_fixture.mel) will overwrite with real cmt output." + ], + "slice": "17A", + "group": "D", + "test": "hybrid_scalar_quat_InterpolationAtTrainingPoints", + "spec": { + "scalar_dim": 1, + "quat_blocks": [ + {"space": "SwingTwist", "axis": [0.0, 1.0, 0.0]} + ] + }, + "N": 3, + "scalar_centers": [[-1.0], [0.0], [1.0]], + "quat_feature_layout": "row = pose, cols = (x, y, z, w)", + "quat_features": [ + [ + [0.0, 0.0, 0.0, 1.0], + [0.0, 0.49999999999999994, 0.0, 0.86602540378443871], + [0.25881904510252074, 0.0, 0.0, 0.96592582628906831] + ] + ], + "axis_angle_closed_forms": { + "pose_0": "identity", + "pose_1": "60 deg around +Y : (0, sin(30 deg), 0, cos(30 deg)) — pure twist about Y", + "pose_2": "30 deg around +X : (sin(15 deg), 0, 0, cos(15 deg)) — swing only (twist identity) about Y" + }, + "targets_primary": [[0.0, 1.0], [0.5, -0.25], [1.0, 0.75]], + "fit_params_primary": { + "kernel": {"type": "Gaussian", "eps": 1.0}, + "poly_degree": -1, + "lambda": 1e-10 + }, + "expectations": { + "status": "OK", + "weights_shape_rows": 9, + "weights_shape_cols_primary": 2, + "feature_norms_size": 1, + "feature_norms_positive": true, + "distance_norm_positive": true, + "residual_norm_upper_bound": 1e-9 + } +} diff --git a/tests/fixtures/cmt_parity_17A_quat_only_full.json b/tests/fixtures/cmt_parity_17A_quat_only_full.json new file mode 100644 index 0000000..b2600e2 --- /dev/null +++ b/tests/fixtures/cmt_parity_17A_quat_only_full.json @@ -0,0 +1,51 @@ +{ + "_comment_disclaimer": [ + "NOT a cmt binary parity fixture.", + "Expected values computed by independent re-implementation per plan section C pseudocode.", + "True cmt binary parity is deferred to Slice 17B (T-30).", + "", + "Slice 17A audit trail: documents the 4-pose Full-mode test input used by", + "tests/test_feature_spec.cpp::HeteroFit_quat_only_FullMode_InterpolatesFixture.", + "Tests hard-code these doubles in C++ to avoid JSON-parse rounding; this file", + "serves as reviewer audit trail and as the payload Slice 17B's cmt binary parity", + "generator (scripts/gen_cmt_fixture.mel) will overwrite with real cmt output." + ], + "slice": "17A", + "group": "C", + "test": "quat_only_FullMode_InterpolatesFixture", + "spec": { + "scalar_dim": 0, + "quat_blocks": [ + {"space": "Full", "axis": [0.0, 0.0, 0.0]} + ] + }, + "N": 4, + "quat_feature_layout": "row = pose, cols = (x, y, z, w)", + "quat_features": [ + [ + [0.0, 0.0, 0.0, 1.0], + [0.0, 0.70710678118654757, 0.0, 0.70710678118654757], + [0.0, 1.0, 0.0, 0.0], + [0.0, 0.38268343236508978, 0.0, 0.92387953251128674] + ] + ], + "axis_angle_closed_forms": { + "pose_0": "identity (0 rad)", + "pose_1": "90 deg around +Y : (0, sin(pi/4), 0, cos(pi/4))", + "pose_2": "180 deg around +Y: (0, 1, 0, 0)", + "pose_3": "45 deg around +Y : (0, sin(pi/8), 0, cos(pi/8))" + }, + "targets": [[0.1], [0.4], [0.9], [0.3]], + "fit_params": { + "kernel": {"type": "Gaussian", "eps": 1.0}, + "poly_degree": -1, + "lambda": 1e-10 + }, + "expectations": { + "status": "OK", + "weights_shape_rows": 4, + "weights_shape_cols": 1, + "residual_norm_upper_bound": 1e-9, + "tolerance_note": "Algorithm-self-consistent tolerance <= 1e-12 per plan Decision 2; residual_norm bound here is the ridge-regression residual gate, not a direct value assertion." + } +} diff --git a/tests/test_feature_spec.cpp b/tests/test_feature_spec.cpp new file mode 100644 index 0000000..b4ec7c2 --- /dev/null +++ b/tests/test_feature_spec.cpp @@ -0,0 +1,635 @@ +// ============================================================================= +// test_feature_spec.cpp — Phase 2A.5 Slice 17A unit tests +// ----------------------------------------------------------------------------- +// 20 TEST blocks across 4 groups: +// A — FeatureSpec / QuatBlock type contract (8 live) +// B — 17A-SCALAR-ORACLE: byte-identical FitResult vs legacy fit() (6 live) +// C — quat-only mixed spaces (Full / SwingTwist) (3 SKIP) +// D — hybrid scalar + quat (3 SKIP) +// +// Groups C and D are the hetero pipeline tests; they land in Slice 17A +// step 3.4 with the full composite distance-matrix builder. They are +// stubbed here with GTEST_SKIP so that the binary compiles in step 3.3 +// but the 17A delivery is not mistaken for complete on Group A + B alone. +// +// Tolerance philosophy (Slice 17A) +// -------------------------------- +// Group B is a bit-identity test — no floating-point tolerance is acceptable. +// std::memcmp over Eigen contiguous storage is the oracle for MatrixX / VectorX +// comparisons. For scalar fields (lambda_used, condition_number, residual_norm, +// kernel.eps, distance_norm), the bit pattern is compared via std::memcpy into +// std::uint64_t, since C++20's std::bit_cast is not available on the Maya 2018 +// GCC 4.8.2 C++11 floor. +// ============================================================================= +#include +#include +#include +#include + +#include +#include + +#include "rbfmax/feature_spec.hpp" +#include "rbfmax/kernel_functions.hpp" +#include "rbfmax/solver.hpp" +#include "rbfmax/types.hpp" + +namespace { + +using rbfmax::FeatureSpec; +using rbfmax::Index; +using rbfmax::KernelParams; +using rbfmax::KernelType; +using rbfmax::MatrixX; +using rbfmax::QuatBlock; +using rbfmax::Scalar; +using rbfmax::SolverSpace; +using rbfmax::Vector3; +using rbfmax::VectorX; +using rbfmax::solver::FitOptions; +using rbfmax::solver::FitResult; +using rbfmax::solver::FitStatus; +using rbfmax::solver::kLambdaAuto; +using rbfmax::solver::SolverPath; + +// --------------------------------------------------------------------------- +// Bit-identity helpers +// --------------------------------------------------------------------------- + +// C++11-compatible replacement for std::bit_cast (C++20). Reads the object +// representation of a Scalar (double) into a 64-bit unsigned integer. +inline std::uint64_t scalar_bits(Scalar x) noexcept { + static_assert(sizeof(Scalar) == sizeof(std::uint64_t), + "Scalar must be 64-bit for bit-identity tests"); + std::uint64_t u = 0; + std::memcpy(&u, &x, sizeof(u)); + return u; +} + +inline bool scalar_bit_equal(Scalar a, Scalar b) noexcept { + return scalar_bits(a) == scalar_bits(b); +} + +// Byte-level equality on Eigen contiguous storage. Checks shape first, then +// memcmp over rows*cols*sizeof(Scalar). Used on MatrixX / VectorX. +bool matrix_bytes_equal(const MatrixX& a, const MatrixX& b) { + if (a.rows() != b.rows()) return false; + if (a.cols() != b.cols()) return false; + if (a.size() == 0) return true; // 0×0 matrices are trivially equal + const std::size_t nbytes = + static_cast(a.size()) * sizeof(Scalar); + return std::memcmp(a.data(), b.data(), nbytes) == 0; +} + +// (vector_bytes_equal removed — 17A oracle has no VectorX field comparison: +// scalars → scalar_bit_equal; matrices → matrix_bytes_equal; +// feature_norms asserts size()==0, never needing byte compare. +// Deletion avoids GCC -Werror=unused-function on Ubuntu CI.) + +// Full 10-field byte-identical comparison of two FitResults on the pre-17A +// surface only. The 4 new tail fields are verified separately — see +// expect_17a_tail_default_with_overlayed_spec. +void expect_fit_result_byte_identical(const FitResult& old_fr, + const FitResult& new_fr) { + // 1. weights — full bit-identity of N × M Eigen storage. + EXPECT_TRUE(matrix_bytes_equal(old_fr.weights, new_fr.weights)) + << "weights byte-identity failed"; + + // 2. poly_coeffs — may be 0×M; the helper handles that. + EXPECT_TRUE(matrix_bytes_equal(old_fr.poly_coeffs, new_fr.poly_coeffs)) + << "poly_coeffs byte-identity failed"; + + // 3. centers — N × D owned copy; must be identical. + EXPECT_TRUE(matrix_bytes_equal(old_fr.centers, new_fr.centers)) + << "centers byte-identity failed"; + + // 4. kernel.type — enum equality. + EXPECT_EQ(static_cast(old_fr.kernel.type), + static_cast(new_fr.kernel.type)); + + // 5. kernel.eps — bit-identity via uint64 pattern, NOT EXPECT_DOUBLE_EQ + // (4-ULP tolerance is a violation of the SCALAR-ORACLE contract). + EXPECT_TRUE(scalar_bit_equal(old_fr.kernel.eps, new_fr.kernel.eps)) + << "kernel.eps bit pattern drift: old=" << scalar_bits(old_fr.kernel.eps) + << " new=" << scalar_bits(new_fr.kernel.eps); + + // 6. poly_degree. + EXPECT_EQ(old_fr.poly_degree, new_fr.poly_degree); + + // 7. lambda_used — bit-identity. + EXPECT_TRUE(scalar_bit_equal(old_fr.lambda_used, new_fr.lambda_used)) + << "lambda_used bit pattern drift"; + + // 8. solver_path enum. + EXPECT_EQ(static_cast(old_fr.solver_path), + static_cast(new_fr.solver_path)); + + // 9. status enum. + EXPECT_EQ(static_cast(old_fr.status), + static_cast(new_fr.status)); + + // 10. condition_number — bit-identity. + EXPECT_TRUE(scalar_bit_equal(old_fr.condition_number, + new_fr.condition_number)) + << "condition_number bit pattern drift"; + + // 11. residual_norm — bit-identity. + EXPECT_TRUE(scalar_bit_equal(old_fr.residual_norm, new_fr.residual_norm)) + << "residual_norm bit pattern drift"; +} + +// Assertion over the 4 Slice 17A tail fields for a scalar-only dispatched +// FitResult. feature_spec must reflect the caller's spec (NOT reset to +// default) — this is the explicit carrier of the heterogeneous API contract, +// even in the scalar-only branch where quat_features is empty. +void expect_17a_tail_default_with_overlayed_spec( + const FitResult& fr, const FeatureSpec& expected_spec) { + // feature_spec: caller's spec is echoed back (overlay in dispatch). + EXPECT_EQ(fr.feature_spec.scalar_dim, expected_spec.scalar_dim); + EXPECT_TRUE(fr.feature_spec.is_scalar_only()); + EXPECT_EQ(fr.feature_spec.quat_blocks.size(), + expected_spec.quat_blocks.size()); + + // quat_features: empty for scalar-only. + EXPECT_TRUE(fr.quat_features.empty()); + + // feature_norms: size 0. + EXPECT_EQ(fr.feature_norms.size(), Index(0)); + + // distance_norm: bit-identical to 0.0. + EXPECT_TRUE(scalar_bit_equal(fr.distance_norm, Scalar(0))) + << "distance_norm must be exactly 0.0 in scalar-only dispatch"; +} + +// --------------------------------------------------------------------------- +// Small deterministic fixtures — no RNG so oracle comparisons stay stable. +// --------------------------------------------------------------------------- + +MatrixX make_centers_2d_4pt() { + MatrixX C(4, 2); + C << 0.0, 0.0, + 1.0, 0.0, + 0.0, 1.0, + 1.0, 1.0; + return C; +} + +MatrixX make_targets_4x1_linear() { + MatrixX Y(4, 1); + Y << 0.0, 1.0, 1.0, 2.0; + return Y; +} + +MatrixX make_targets_4x3_multi() { + MatrixX Y(4, 3); + Y << 0.0, 0.5, -0.25, + 1.0, 0.5, 0.75, + 1.0, -0.5, 0.75, + 2.0, -0.5, 1.75; + return Y; +} + +} // namespace + +// ============================================================================ +// Group A — FeatureSpec / QuatBlock type contract +// ============================================================================ + +TEST(FeatureSpec, DefaultIsScalarOnly) { + FeatureSpec spec; + EXPECT_EQ(spec.scalar_dim, Index(0)); + EXPECT_TRUE(spec.quat_blocks.empty()); + EXPECT_TRUE(spec.is_scalar_only()); +} + +TEST(FeatureSpec, ScalarOnlyFromDimCtor) { + FeatureSpec spec(Index(3)); + EXPECT_EQ(spec.scalar_dim, Index(3)); + EXPECT_TRUE(spec.is_scalar_only()); + EXPECT_EQ(spec.total_distance_columns(Index(5)), Index(5)); // N cols +} + +TEST(FeatureSpec, TotalDistanceColumns_ScalarOnly) { + // scalar_dim > 0 => N columns; scalar_dim == 0 => 0 columns. + FeatureSpec with_scalars(Index(4)); + EXPECT_EQ(with_scalars.total_distance_columns(Index(7)), Index(7)); + + FeatureSpec without; + EXPECT_EQ(without.total_distance_columns(Index(7)), Index(0)); +} + +TEST(FeatureSpec, TotalDistanceColumns_MixedSpaces) { + // Full → 1 col per pose; Swing/Twist/SwingTwist → 2 cols per pose. + // Spec: scalar_dim=2, 3 quat blocks (Full, Swing, SwingTwist). + // For N=5: 5 (scalar) + 1*5 (Full) + 2*5 (Swing) + 2*5 (SwingTwist) = 30. + FeatureSpec spec(Index(2)); + spec.quat_blocks.push_back(QuatBlock(SolverSpace::Full, Vector3::Zero())); + spec.quat_blocks.push_back( + QuatBlock(SolverSpace::Swing, Vector3::UnitY())); + spec.quat_blocks.push_back( + QuatBlock(SolverSpace::SwingTwist, Vector3::UnitY())); + + EXPECT_EQ(spec.total_distance_columns(Index(5)), Index(5 + 5 + 10 + 10)); + EXPECT_FALSE(spec.is_scalar_only()); +} + +TEST(QuatBlock, DefaultAxisIsZero) { + QuatBlock qb; + EXPECT_EQ(qb.space, SolverSpace::Full); + EXPECT_TRUE(qb.axis.isZero(Scalar(0))); // exact zero +} + +TEST(QuatBlock, SolverSpaceEnumRoundTrip) { + // Verifies the enum backing type and the four canonical values. + EXPECT_EQ(static_cast(SolverSpace::Full), 0); + EXPECT_EQ(static_cast(SolverSpace::Swing), 1); + EXPECT_EQ(static_cast(SolverSpace::Twist), 2); + EXPECT_EQ(static_cast(SolverSpace::SwingTwist), 3); + + EXPECT_EQ(FeatureSpec::cols_per_pose(SolverSpace::Full), Index(1)); + EXPECT_EQ(FeatureSpec::cols_per_pose(SolverSpace::Swing), Index(2)); + EXPECT_EQ(FeatureSpec::cols_per_pose(SolverSpace::Twist), Index(2)); + EXPECT_EQ(FeatureSpec::cols_per_pose(SolverSpace::SwingTwist), Index(2)); +} + +TEST(FeatureSpec, NonEmptyQuatBlocksImpliesNonScalarOnly) { + FeatureSpec spec(Index(0)); + spec.quat_blocks.push_back(QuatBlock(SolverSpace::Full, Vector3::Zero())); + EXPECT_FALSE(spec.is_scalar_only()); +} + +TEST(FeatureSpec, MoveConstructible) { + // std::vector inside FitResult requires FeatureSpec to survive + // container operations. C++11 nothrow-moveable wins the "is_nothrow_move + // _constructible" trait check when QuatBlock is trivially relocatable. + FeatureSpec src(Index(4)); + src.quat_blocks.push_back( + QuatBlock(SolverSpace::SwingTwist, Vector3::UnitZ())); + + FeatureSpec moved(std::move(src)); + EXPECT_EQ(moved.scalar_dim, Index(4)); + EXPECT_EQ(moved.quat_blocks.size(), std::size_t(1)); + EXPECT_EQ(moved.quat_blocks[0].space, SolverSpace::SwingTwist); +} + +// ============================================================================ +// Group B — 17A-SCALAR-ORACLE (hard gate) +// +// For each of 6 (kernel × λ × poly) combinations, fit both via legacy +// solver::fit(C, Y, opts, λ) and the new heterogeneous solver::fit(C, {}, Y, +// opts, FeatureSpec(C.cols()), λ). Assert: +// * all 10 pre-17A FitResult fields are byte-identical +// * the 4 new tail fields are default-constructed except feature_spec, +// which must reflect the caller's (scalar_dim=C.cols(), empty blocks). +// ============================================================================ + +TEST(Oracle17A, scalar_only_baseline_ByteIdentical_Gaussian) { + const MatrixX C = make_centers_2d_4pt(); + const MatrixX Y = make_targets_4x1_linear(); + FitOptions opt(KernelParams(KernelType::kGaussian, 1.0)); + opt.poly_degree = -1; + + const FitResult old_fr = + rbfmax::solver::fit(C, Y, opt, Scalar(1e-8)); + const FeatureSpec spec(C.cols()); + const FitResult new_fr = rbfmax::solver::fit( + C, std::vector{}, Y, opt, spec, Scalar(1e-8)); + + expect_fit_result_byte_identical(old_fr, new_fr); + expect_17a_tail_default_with_overlayed_spec(new_fr, spec); + EXPECT_EQ(new_fr.status, FitStatus::OK); +} + +TEST(Oracle17A, scalar_only_baseline_ByteIdentical_ThinPlateSpline) { + const MatrixX C = make_centers_2d_4pt(); + const MatrixX Y = make_targets_4x3_multi(); + FitOptions opt(KernelParams(KernelType::kThinPlateSpline, 1.0)); + opt.poly_degree = 1; + + const FitResult old_fr = + rbfmax::solver::fit(C, Y, opt, Scalar(1e-10)); + const FeatureSpec spec(C.cols()); + const FitResult new_fr = rbfmax::solver::fit( + C, std::vector{}, Y, opt, spec, Scalar(1e-10)); + + expect_fit_result_byte_identical(old_fr, new_fr); + expect_17a_tail_default_with_overlayed_spec(new_fr, spec); +} + +TEST(Oracle17A, scalar_only_baseline_ByteIdentical_Cubic_WithPoly) { + const MatrixX C = make_centers_2d_4pt(); + const MatrixX Y = make_targets_4x1_linear(); + FitOptions opt(KernelParams(KernelType::kCubic, 1.0)); + opt.poly_degree = 1; + + const FitResult old_fr = + rbfmax::solver::fit(C, Y, opt, Scalar(1e-10)); + const FeatureSpec spec(C.cols()); + const FitResult new_fr = rbfmax::solver::fit( + C, std::vector{}, Y, opt, spec, Scalar(1e-10)); + + expect_fit_result_byte_identical(old_fr, new_fr); + expect_17a_tail_default_with_overlayed_spec(new_fr, spec); + EXPECT_GT(new_fr.poly_coeffs.rows(), 0); // poly tail actually present +} + +TEST(Oracle17A, scalar_only_baseline_ByteIdentical_LambdaAuto) { + const MatrixX C = make_centers_2d_4pt(); + const MatrixX Y = make_targets_4x1_linear(); + FitOptions opt(KernelParams(KernelType::kGaussian, 1.0)); + opt.poly_degree = -1; + + const FitResult old_fr = rbfmax::solver::fit(C, Y, opt, kLambdaAuto); + const FeatureSpec spec(C.cols()); + const FitResult new_fr = rbfmax::solver::fit( + C, std::vector{}, Y, opt, spec, kLambdaAuto); + + expect_fit_result_byte_identical(old_fr, new_fr); + expect_17a_tail_default_with_overlayed_spec(new_fr, spec); + // GCV should have selected a positive λ. + EXPECT_GT(new_fr.lambda_used, Scalar(0)); +} + +TEST(Oracle17A, scalar_only_baseline_PreservesFitStatus) { + // Trigger INSUFFICIENT_SAMPLES via an empty-rows center matrix; the + // dispatch must propagate legacy's validate_inputs status verbatim. + const MatrixX C(0, 2); // 0 samples + const MatrixX Y(0, 1); + FitOptions opt(KernelParams(KernelType::kGaussian, 1.0)); + + const FitResult old_fr = + rbfmax::solver::fit(C, Y, opt, Scalar(1e-8)); + const FeatureSpec spec(C.cols()); + const FitResult new_fr = rbfmax::solver::fit( + C, std::vector{}, Y, opt, spec, Scalar(1e-8)); + + EXPECT_EQ(old_fr.status, FitStatus::INSUFFICIENT_SAMPLES); + EXPECT_EQ(new_fr.status, FitStatus::INSUFFICIENT_SAMPLES); + EXPECT_EQ(old_fr.solver_path, SolverPath::FAILED); + EXPECT_EQ(new_fr.solver_path, SolverPath::FAILED); + // Spec is still overlaid even on the INSUFFICIENT_SAMPLES path. + expect_17a_tail_default_with_overlayed_spec(new_fr, spec); +} + +TEST(Oracle17A, scalar_only_baseline_PreservesSolverPath) { + // A well-conditioned small system normally takes the LLT path. Confirm + // dispatch preserves solver_path verbatim across a mix-output case. + const MatrixX C = make_centers_2d_4pt(); + const MatrixX Y = make_targets_4x3_multi(); + FitOptions opt(KernelParams(KernelType::kGaussian, 1.0)); + opt.poly_degree = -1; + + const FitResult old_fr = + rbfmax::solver::fit(C, Y, opt, Scalar(1e-6)); + const FeatureSpec spec(C.cols()); + const FitResult new_fr = rbfmax::solver::fit( + C, std::vector{}, Y, opt, spec, Scalar(1e-6)); + + expect_fit_result_byte_identical(old_fr, new_fr); + expect_17a_tail_default_with_overlayed_spec(new_fr, spec); + // Weight shape sanity: N=4 centers × M=3 outputs. + EXPECT_EQ(new_fr.weights.rows(), Index(4)); + EXPECT_EQ(new_fr.weights.cols(), Index(3)); +} + +// --------------------------------------------------------------------------- +// Group C/D fixtures (Slice 17A step 3.4) +// +// The two JSON files under tests/fixtures/ document these hand-computed +// quaternion inputs verbatim; tests hard-code the doubles here for +// numerical stability (avoids JSON parse rounding). The fixtures serve as +// audit trail and as the JSON payload 17B's cmt-binary-parity generator +// will overwrite. +// +// Quat storage convention in MatrixX: row = pose, columns = (x, y, z, w). +// +// Wrapped in an anonymous namespace to keep internal linkage — these +// helpers are translation-unit-local and must not leak to `rbfmax::` +// (GCC -Werror=missing-declarations catches bare file-scope functions +// without a prior declaration; MSVC does not). +// --------------------------------------------------------------------------- + +namespace { + +// Full-mode fixture: N=4, axis-angle closed-form unit quaternions. +// pose 0 : identity → (0, 0, 0, 1) +// pose 1 : 90° around +Y → (0, sin(π/4), 0, cos(π/4)) = (0, √2/2, 0, √2/2) +// pose 2 : 180° around +Y → (0, 1, 0, 0) +// pose 3 : 45° around +Y → (0, sin(π/8), 0, cos(π/8)) +MatrixX make_quat_fixture_full_N4() { + MatrixX Q(4, 4); + const Scalar sqrt2_over_2 = Scalar(0.70710678118654757); // √2/2 + const Scalar sin_pi_8 = Scalar(0.38268343236508978); // sin(π/8) + const Scalar cos_pi_8 = Scalar(0.92387953251128674); // cos(π/8) + Q << 0.0, 0.0, 0.0, 1.0, // pose 0 + 0.0, sqrt2_over_2, 0.0, sqrt2_over_2, // pose 1 + 0.0, 1.0, 0.0, 0.0, // pose 2 + 0.0, sin_pi_8, 0.0, cos_pi_8; // pose 3 + return Q; +} + +// SwingTwist-mode fixture (axis=Y): N=3, unit quaternions with non-trivial +// swing and twist components. +// pose 0 : identity +// pose 1 : 60° around +Y (pure twist around Y) +// pose 2 : 30° around +X (pure swing away from Y axis) +MatrixX make_quat_fixture_swingtwist_N3() { + MatrixX Q(3, 4); + const Scalar s30 = Scalar(0.25881904510252074); // sin(15°) + const Scalar c30 = Scalar(0.96592582628906831); // cos(15°) + const Scalar s60 = Scalar(0.49999999999999994); // sin(30°) + const Scalar c60 = Scalar(0.86602540378443871); // cos(30°) + Q << 0.0, 0.0, 0.0, 1.0, // identity + 0.0, s60, 0.0, c60, // 60° about Y (twist only for axis=Y) + s30, 0.0, 0.0, c30; // 30° about X (swing-dominant for axis=Y) + return Q; +} + +} // namespace + +// ============================================================================ +// Group C — quat-only mixed spaces (Slice 17A step 3.4) +// ============================================================================ + +TEST(HeteroFit, quat_only_FullMode_InterpolatesFixture) { + // Heterogeneous fit with one Full-mode quat block (N=4, B=1, scalar_dim=0). + // Expected: status OK, weights shape = (total_cols × M) = (4 × 1), + // residual_norm small (ridge with λ=1e-8 on well-conditioned 4-pose data). + const MatrixX scalar_centers(0, 0); + std::vector quat_features; + quat_features.push_back(make_quat_fixture_full_N4()); + MatrixX Y(4, 1); + Y << 0.1, 0.4, 0.9, 0.3; + + FeatureSpec spec; + spec.quat_blocks.push_back(QuatBlock(SolverSpace::Full, Vector3::Zero())); + + FitOptions opt(KernelParams(KernelType::kGaussian, 1.0)); + opt.poly_degree = -1; + + // λ=1e-10: for ridge regression, residual_norm scales as O(λ/(σ_min+λ)), + // so to hit the 1e-9 interpolation gate we need λ < 1e-9. Matches the + // hybrid Group D calibration. + const FitResult fr = rbfmax::solver::fit( + scalar_centers, quat_features, Y, opt, spec, Scalar(1e-10)); + + ASSERT_EQ(fr.status, FitStatus::OK) + << "hetero Full-mode fit must succeed on the 4-pose fixture"; + EXPECT_NE(fr.solver_path, SolverPath::FAILED); + + // Full-mode column count: 1 per pose × N = 4 cols total. + EXPECT_EQ(fr.weights.rows(), Index(4)) + << "Full-mode weights.rows() == total_distance_columns(N) == N"; + EXPECT_EQ(fr.weights.cols(), Index(1)); + + // Tail fields populated (distinguishes hetero branch from scalar-only). + EXPECT_FALSE(fr.feature_spec.is_scalar_only()); + EXPECT_EQ(fr.quat_features.size(), std::size_t(1)); + EXPECT_EQ(fr.quat_features[0].rows(), Index(4)); + EXPECT_EQ(fr.quat_features[0].cols(), Index(4)); + + // Interpolation property via ridge residual (tol 1e-9 per plan). + EXPECT_LT(fr.residual_norm, Scalar(1e-9)) + << "ridge-regression residual too large; check kernel / lambda"; +} + +TEST(HeteroFit, quat_only_SwingTwistMode_CorrectColumnCount) { + // Column count invariant: SwingTwist contributes 2N per block. + // N=3, B=1 → total_cols = 2*3 = 6. + const MatrixX scalar_centers(0, 0); + std::vector quat_features; + quat_features.push_back(make_quat_fixture_swingtwist_N3()); + MatrixX Y(3, 1); + Y << 1.0, -2.0, 3.0; + + FeatureSpec spec; + spec.quat_blocks.push_back( + QuatBlock(SolverSpace::SwingTwist, Vector3::UnitY())); + + // Helper's own arithmetic verified first — avoids surprise if the fit + // pipeline silently renormalises total_distance_columns(). + EXPECT_EQ(spec.total_distance_columns(Index(3)), Index(6)); + + FitOptions opt(KernelParams(KernelType::kGaussian, 1.0)); + const FitResult fr = rbfmax::solver::fit( + scalar_centers, quat_features, Y, opt, spec, Scalar(1e-8)); + + ASSERT_EQ(fr.status, FitStatus::OK); + EXPECT_EQ(fr.weights.rows(), Index(6)) + << "SwingTwist weights.rows() == 2 * N"; + EXPECT_EQ(fr.weights.cols(), Index(1)); +} + +TEST(HeteroFit, quat_only_ValidationRejectsBlockMismatch) { + // spec declares 2 blocks but caller supplies only 1 quat_features entry. + const MatrixX scalar_centers(0, 0); + std::vector quat_features; + quat_features.push_back(make_quat_fixture_full_N4()); + MatrixX Y(4, 1); + Y << 0.1, 0.2, 0.3, 0.4; + + FeatureSpec spec; + spec.quat_blocks.push_back(QuatBlock(SolverSpace::Full, Vector3::Zero())); + spec.quat_blocks.push_back( + QuatBlock(SolverSpace::SwingTwist, Vector3::UnitY())); // extra block + + FitOptions opt(KernelParams(KernelType::kGaussian, 1.0)); + const FitResult fr = rbfmax::solver::fit( + scalar_centers, quat_features, Y, opt, spec, Scalar(1e-8)); + + EXPECT_EQ(fr.status, FitStatus::INVALID_INPUT) + << "block count mismatch must surface INVALID_INPUT"; + EXPECT_EQ(fr.solver_path, SolverPath::FAILED); +} + +// ============================================================================ +// Group D — hybrid scalar + quat (Slice 17A step 3.4) +// ============================================================================ + +TEST(HeteroFit, hybrid_scalar_quat_FitStatusOK_OnValidInput) { + // N=3, scalar_dim=1, B=1 SwingTwist → cols = 3 + 2*3 = 9. + MatrixX scalar_centers(3, 1); + scalar_centers << -1.0, 0.0, 1.0; + std::vector quat_features; + quat_features.push_back(make_quat_fixture_swingtwist_N3()); + MatrixX Y(3, 1); + Y << 0.5, 1.5, 2.5; + + FeatureSpec spec(Index(1)); + spec.quat_blocks.push_back( + QuatBlock(SolverSpace::SwingTwist, Vector3::UnitY())); + + FitOptions opt(KernelParams(KernelType::kGaussian, 1.0)); + const FitResult fr = rbfmax::solver::fit( + scalar_centers, quat_features, Y, opt, spec, Scalar(1e-8)); + + ASSERT_EQ(fr.status, FitStatus::OK); + EXPECT_EQ(fr.weights.rows(), Index(9)); + EXPECT_EQ(fr.weights.cols(), Index(1)); + + // Feature-norms captured for the single scalar column. + EXPECT_EQ(fr.feature_norms.size(), Index(1)); + EXPECT_GT(fr.feature_norms(0), Scalar(0)); + // Distance norm (scalar-block Frobenius) populated non-zero. + EXPECT_GT(fr.distance_norm, Scalar(0)); +} + +TEST(HeteroFit, hybrid_scalar_quat_InterpolationAtTrainingPoints) { + // Same hybrid fit as above, residual gate at 1e-9. For ridge regression + // with very small λ this approximates the minimum-norm solution; on a + // well-conditioned N=3 system the residual is dominated by λ, not + // conditioning, so the tolerance is comfortable. + MatrixX scalar_centers(3, 1); + scalar_centers << -1.0, 0.0, 1.0; + std::vector quat_features; + quat_features.push_back(make_quat_fixture_swingtwist_N3()); + MatrixX Y(3, 2); + Y << 0.0, 1.0, + 0.5, -0.25, + 1.0, 0.75; + + FeatureSpec spec(Index(1)); + spec.quat_blocks.push_back( + QuatBlock(SolverSpace::SwingTwist, Vector3::UnitY())); + + FitOptions opt(KernelParams(KernelType::kGaussian, 1.0)); + const FitResult fr = rbfmax::solver::fit( + scalar_centers, quat_features, Y, opt, spec, Scalar(1e-10)); + + ASSERT_EQ(fr.status, FitStatus::OK); + EXPECT_LT(fr.residual_norm, Scalar(1e-9)) + << "hybrid ridge residual exceeds 1e-9 at training points"; +} + +// Compile-time verification that Slice 17A's FitResult tail additions have +// NOT broken the noexcept lifecycle of FitResult itself — the invariant that +// propagates into solver::predict's noexcept contract. Constraint #6 from +// the Step 3.4 plan asks us to catch a silent regression where (e.g.) adding +// std::vector causes FitResult default ctor / move / destructor to +// silently drop noexcept; the trait checks here are the tightest portable +// signal under C++11 (the expression-form noexcept(solver::predict(fr, x)) +// depends on MSVC's Eigen::Ref construction being noexcept, which is not +// guaranteed across compiler / Eigen-version combinations and therefore +// unusable as a hard gate). The predict functions themselves carry the +// noexcept keyword in solver.hpp and that remains the canonical guarantee. +TEST(HeteroFit, hybrid_scalar_quat_PredictNoexcept) { + static_assert( + std::is_nothrow_default_constructible::value, + "FitResult default ctor must remain noexcept across Slice 17A tail fields"); + static_assert( + std::is_nothrow_destructible::value, + "FitResult destructor must remain noexcept across Slice 17A tail fields"); + static_assert( + std::is_nothrow_default_constructible::value, + "FeatureSpec default ctor must be noexcept (17A contract)"); + static_assert( + std::is_nothrow_destructible::value, + "FeatureSpec destructor must be noexcept"); + static_assert( + std::is_nothrow_default_constructible::value, + "QuatBlock default ctor must be noexcept"); + // The function-level noexcept keyword on solver::predict / predict_scalar + // in solver.hpp remains the canonical guarantee; this TEST block validates + // the FitResult-side preconditions that the keyword relies on. + SUCCEED(); +}