CultLib/docs/cultlib-gaps-map.md
Meta Vens 50512e6e4f cultlib-gaps map: BLS early-stop diagnosis and CultMath rotation algebra
Imagination pass (session self-2026-10-10-ag): the cost-progress stop at
BoundedLeastSquares.cs:275 reports Converged with a bound still violated;
fix follows Lawson-Hanson. Quaternion algebra surface and the
mul(float4,float4) overload trap.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-10 07:12:09 +02:00

18 KiB

CultLib Gaps: Map

Campaign cultlib-gaps (target: docs/cultlib-gaps-target.md). This file holds body facts, probes and rationale for the cuts the mind's cut_specs name; the specs cite it by section and fact number.

Pinned base: origin/main at 6bdc23c9 (2026-10-10). Imagination pass imagination-cl-gaps-cultmath, session self-2026-10-10-ag. Probes live under the session scratchpad, imagination-cl-gaps-cultmath/ (probe1 traced copy of the solver and the warm/cold sweep, suite a scratch copy of packages/cultmath with the candidate fix, probe2 the sweep against that copy, probe3 the overload probe). None of it is committed.

Cut order

  1. bls-release-stop: BoundedLeastSquares stops on a release only when the release is not real.
  2. bls-input-field: a validation helper names the invalid input (follow-up bls-invalid-input-field). Blocked on question cultlib-gaps:question:bls-input-error-shape.
  3. cultmath-release-0-4-1: release org.gamecult.cultmath with cuts 1 and 2. Version follows question cultlib-gaps:question:bls-fix-release-lane.
  4. cultmath-rotation-algebra: quaternion algebra, swing-twist and joint clamps (ruling aetheria-release:ruling:solver-in-gamecult-animation). Ships in the CultMath release that the first GameCult.Animation release pins; that release is mapped with GameCult.Animation.

B. BoundedLeastSquares stopping

Finding aetheria-release:finding:cut-thrust-allocator-core.s1.solver-converged-early (Soul, Medium): a cold start reports Converged at a non-optimal point, so Aetheria's thrust allocator (tolerance 1e-12) answers differently warm and cold, by up to .078 throttle.

Body facts

  • B1. One runtime ships the solver. BoundedLeastSquares exists only in packages/cultmath/src/CultMath/BoundedLeastSquares.cs (C#, netstandard2.1 and net8.0). The type doc says there is no HLSL mirror; grep over packages/* finds no Rust, TypeScript, Python or Kotlin sibling, and no other CultLib project calls it. Its consumer outside CultLib is Aetheria's ThrustAllocator.
  • B2. The defect reproduces against CultMath at base, outside Aetheria. probe1 rebuilds Soul's allocator matrix (Soul's ProbeAllocator2, constants H=100, Wy=1000, W=.3, kappa=10, rho=.1, EMax=2) for hull Synthetic(8, 7005) and the 8th intent of System.Random(5), (-0.8505, 0.6739, 0.5587): 15 rows, 12 columns. Real BoundedLeastSquares.Solve at tolerance 1e-12 from x = 0: Converged, 13 iterations, cost 36.0313737, x7 = 1. Re-solving from that answer: 1 iteration, cost 36.0312608, x7 = 0.9222145, which is also the warm answer. The same literals (dumped with "R" floats, probe1 dump) reproduce it at tolerance 1e-9 and 1e-12, and not at the default 1e-6, where the cold solve stops at iteration 10 inside the tolerance and the re-solve takes 0 iterations.
  • B2a. The W5 problem as literals (probe1 dump, float "R" round trip; m = 15, n = 12, row-major; the regression fixture for bls-release-stop). Nonzero entries of A:
    • row 0: a[0,2] = 100f, a[0,8] = -100f, a[0,9] = 100f
    • row 1: 13.882518f, 20.86695f, 0, 14.8965845f, 12.38046f, 22.115599f, -13.706654f, 15.85789f, 0, 0, -100f, 100f
    • row 2: -187.72859f, 239.31168f, 158.87828f, -270.31165f, 235.38467f, -14.704258f, 121.61232f, 244.81305f, 0, 0, 0, 0
    • rows 3-6: a[3+k, 8+k] = 0.3f (k = 0..3); rows 7-14: a[7+i, i] = 0.1f (i = 0..7)
    • b = {0, 67.385445f, 558.70917f, -3f, -3f, -3f, -3f, 0 x 8}; lo = 0; hi = 1f for columns 0-7, 2f for 8-11. At base, tolerance 1e-12 from x = 0: Converged, 13 iterations, x[7] = 1, cost 36.031373704487095; a re-solve from that answer: 1 iteration, x[7] = 0.92221445f, cost 36.031260771352414.
  • B3. The stopping rule that fires is the progress rule, after a release (BoundedLeastSquares.cs:273-278). A traced copy of the solver, bit-equal to the real one on this problem, shows the last two passes:
    • iteration 12: full step lands; column 2 sits at its lower bound with gradient -0.00205 (violation 2.05e-3, tolerance 1.5e-7); the solver releases it and records fMark = f.
    • iteration 13: the re-solve moves x2 from 0 to 2.05e-7. The cost falls by about 2.1e-10, and the stall band 1e-15 * fScale is 2.11e-10, so f >= fMark - band holds and the solver reports Converged. Column 7 is then at its upper bound with gradient +0.00145, ten thousand times the tolerance; releasing it is worth 1.1e-4 of cost.
  • B4. Why the progress rule cannot judge a release. Near a minimiser, a cost difference resolves only about the square root of what a gradient comparison resolves: a real step of size dx against gradient g lowers the cost by about gdx, while the cost itself is a sum of terms of size fScale. Here gdx (2e-10) sat exactly at the cost's rounding band. The gradient comparison is far from its own rounding: for column 7, the bound on the double gradient's rounding error, (n+1) * eps * (|A^T b|_7 + sum_k |(A^T A)_7k x_k|) (Higham's inner-product bound, gamma_(n+1)), is 9.3e-10, and the violation is 1.56e6 times that.
  • B5. What the progress rule is for (commit 43035bc8, and the tests at BoundedLeastSquaresTests.cs:774-942). It stops two kinds of pass that can never reach the tolerance: a free-set re-solve whose residual is rounding on a near-singular free set, and a release whose violation is itself rounding (kktTol at or below the noise, including kktTol = 0 when the box point nearest the origin is stationary). Deleting it outright breaks both: with a noise gate on every test and no progress rule, 8 of 66 BLS tests fail (wide, integer-tie, collinear, cond spin to the cap; NearSingularFreeSetsConfirmOnceAndStop sees 0 confirmations). Exempting releases from the progress rule with no noise gate fails one: StalledPassesConvergeInsteadOfSpinningToTheCap("plain") seed 141, b = 0, kktTol = 0, cycles release 3 (violation 1.7e-18) / stationarity miss / block on 3 until the cap.
  • B6. The candidate fix, measured. In a scratch copy of the package (suite, diff suite/clean.diff, 44 lines): (a) a bound column is a release candidate only if its violation exceeds both kktTol and the rounding bound of B4; (b) the progress rule judges only passes that release nothing; a pass that releases sets fMark to +Inf, so a release is never measured against a cost recorded before it. Results:
    • all 66 tests in BoundedLeastSquaresTests pass, including DefaultSolvesAreBitIdenticalToTheirGoldenBits and every stall and noise-floor test; the same holds with the bound's factor at 1, 16 and 256 in place of n+1;
    • the W5 problem: cold Converged in 14 iterations at the warm answer (x7 = 0.92221, cost 36.03126077);
    • 400 synthetic hulls x 60 random intents (24,000 calls, warm and cold) at tolerance 1e-12: every call Converged, max iterations 45, warm and cold agree within 6e-8 on every call (the original: one call off by .078 in the first 2,400). At the default tolerance warm and cold still differ, by more than .1 in 10,273 of 24,000 calls: that is the documented resolution of 1e-6 times the gradient scale, pinned by AnAllocatorShapedProblemIsStartIndependentAtATightTolerance (spread >= 0.5), and not this cut's subject.
  • B7. Released columns the factor drops. With the fix, 266 releases in the suite land on a column that SolveFreeSet marks dependent (p = 0). The next pass is then a stationarity miss with no release, and the progress rule stops it, as before. No release in the suite or the 24,000 allocator calls produced a step that moves the released column outward. A Lawson-Hanson style rejection guard (re-bind the column, try the next candidate) was prototyped and is not needed by any measured case; it is not in the cut.

Prior art, and which one the solver follows

  • The solver is a primal active-set method on the normal equations: solve the free set, step to the first bound crossed, release the most violated bound only after a full step lands. That is the outer/inner loop of Lawson-Hanson NNLS (1974) generalised to two-sided bounds, which is the shape of Stark-Parker BVLS (1995) and of scipy's lsq_linear(method='bvls'). It differs in factorisation (Cholesky of A^T A, not QR of A) and in stopping.
  • Lawson-Hanson stop on the dual test only (netlib NNLS, label 350: IF (WMAX .le. ZERO)), and guard roundoff per candidate: a column whose new norm is negligible against rounding (DIFF(UNORM+...)) or whose new coefficient is not positive (ZTEST .gt. ZERO) is rejected with W(J)=ZERO and the next candidate tried. There is no cost-based stop.
  • scipy's BVLS (scipy/optimize/_lsq/bvls.py, main) measures KKT exactly as CultMath's status doc states (g * on_bound for bound variables, |g| for free ones, max over components) and has two stops: optimality below tol (status 1) and a relative cost change below tol (status 2), reported as different statuses.
  • CultMath reported its cost stop as Converged, which its status doc defines as "the KKT conditions hold". The fix follows Lawson-Hanson: a release is judged on the gradient against its rounding (the role of their DIFF test), and the cost stop is kept only where no gradient test can decide (a free-set residual at rounding level), which is where scipy's status 2 also lives.

Rejected alternatives

  • Remove the progress rule and stop on the gradient alone: B5, 8 failures.
  • Report the cost stop as a new status (scipy's status 2): the W5 answer would still be wrong, just labelled, and a new enum member is a breaking change for callers that switch on the status.
  • Measure progress from the step (the exact model decrease t(1-t/2)|L^-1 g|^2) instead of differencing the cost: exact, but a noise-level gradient still yields a tiny positive decrease, so it needs the same noise threshold, on a quantity further from the KKT test than the gradient itself.

I. Naming the invalid input

Follow-up cultlib-gaps:follow_up:bls-invalid-input-field (from Soul finding cut-bls-kkt-tolerance.s1.invalid-input-names-no-field): every validation failure returns the same InvalidInput.

  • I1. Validation today is three sites in Solve: tolerance (:147-148), bounds (:149-151), and finiteness of A and b folded into the A^T b build (:153-161), which cannot tell A from b.
  • I2. Precedent: LAPACK reports a bad argument as INFO = -i, naming the argument and never its value. .NET throws ArgumentException with a parameter name, which the solver's no-throw, zero-allocation contract rules out.
  • Recommended shape (question bls-input-error-shape, option validate-helper): a public BoundedLeastSquares.Validate(m, n, a, b, lo, hi, kktRelativeTolerance) returning a new enum BoundedLeastSquaresInput { Valid, InvalidTolerance, InvalidBounds, NonFiniteMatrix, NonFiniteTarget }, the first failure in that order. Solve calls it and returns InvalidInput for anything but Valid, so validation has one owner and the folded finiteness check in the A^T b loop is deleted. Additive: InvalidInput and every existing result stay.

R. Release 0.4.1

  • R1. org.gamecult.cultmath is at 0.4.0 (package.json:4, CultMath.csproj:13, tag cultmath-unity-v0.4.0). scripts/check-changelog-semver.mjs:97-101: in 0.y.z a minor bump is the breaking lane; additions and fixes take the patch lane (0.2.4 added erf as a patch).
  • R2. Cut 1 changes results only where the old solver stopped after a release whose violation was above both the caller's tolerance and the gradient's rounding, that is, where it reported Converged in breach of its status doc. The default path's golden bits are unchanged (B6). The 0.4.0 precedent (ruling cultmath-0-4-0-breaking-r2) declared the snoise tie fix Breaking because two valid functions gave different values; here the old value broke the documented contract. Question bls-fix-release-lane asks which lane; the recommendation is fixed-patch (0.4.1, ### Fixed).
  • The release procedure is the 0.4.0 cut's (cultlib-gaps:cut_spec:cut-cultmath-release.r1): Windows-only DLL rebuild on Starfire, tests on Yggdrasil, public-surface diff, Unity import smoke, Soul before the tag.

Q. Rotation algebra

Ruling aetheria-release:ruling:solver-in-gamecult-animation: CultMath gains quaternion algebra, swing-twist and joint clamps; the aim-chain solve lives in a new package GameCult.Animation. The interface the consumer needs is the Aetheria map's section "The rig's interface to GameCult.Animation and CultMath" (Aetheria origin/master 30d02719, docs/aetheria-release-map.md:7559-7601).

Body facts

  • Q1. What exists. quaternion (src/CultMath/quaternion.cs, 95 lines): public float fields x, y, z, w; identity; LookRotation(float3 forward, float3 up); implicit conversions to and from float4; scalar ==. math.normalize(quaternion) (math.cs:225-235) returns identity below length 1e-20. No mul, rotate, conjugate, inverse, axis-angle, slerp, swing-twist or clamp exists. math is a partial class (math.cs, math.Interval.cs, math.Affine.cs, math.Phacelle.cs); math.rotate(float2, float) exists (math.cs:240).
  • Q2. Serialization. GameCult.Caching.MessagePack serialises quaternion as the float array [x, y, z, w] (CultMathResolver.cs:170-171) and as an xyzw object in JSON (CultMathJson.cs:468). No other runtime (cultcache-rs, -ts, -py) has a quaternion type. The cut adds static functions only: no field, layout or formatter changes, so no wire change and no cross-runtime fixture.
  • Q3. Runtimes. The algebra is C#-only, like LookRotation: shaders carry rotations as float4 (the type's doc), the mirror tests map shader functions to C# and never the reverse (HlslSourceCompatibilityTests.cs:114-117), and the native core (native/cultmath-core) is a Voronoi tone helper. The Unity bridge already converts quaternion (UnityConversions.cs:16,23); it needs no change.
  • Q4. An overload trap, probed. quaternion converts implicitly from float4, and CultMath has no mul(float4, float4). probe3 declares mul(quaternion, quaternion) beside scalar and matrix mul overloads and calls it with two float4: it compiles and returns the quaternion product. HLSL's mul(vector, vector) is the inner product, and Unity.Mathematics 1.3.3 (the copy in Aetheria's PackageCache) defines mul(float4 a, float4 b) as such (matrix.gen.cs:129), so the same call there is a dot product. The cut therefore adds HLSL's vector mul(float2|3|4, float2|3|4) (dot) beside the quaternion mul, so a float4 pair keeps its shader meaning.
  • Q5. Unity.Mathematics as reference, not template. Its signatures for mul, rotate, conjugate, inverse, AxisAngle and slerp (quaternion.cs:99,493,502,630,651,678) are the consumer's stated convention ("mul(a, b) applies b first"). CultMath's design rule is to state semantics in code and not copy Unity.Mathematics (docs/design.md:465-469), so the formulas are written from the definitions below.

Semantics the cut lands

Quaternions are Hamilton, stored xyzw, unit for rotation; angles are radians, right-handed about the axis.

  • math.mul(quaternion a, quaternion b): the Hamilton product; rotate(mul(a, b), v) == rotate(a, rotate(b, v)).
  • math.rotate(quaternion q, float3 v): t = 2 cross(q.xyz, v); v + q.w t + cross(q.xyz, t); q assumed unit.
  • math.conjugate(q): (-x, -y, -z, w). math.inverse(q): conjugate(q) / dot(q, q), no guard (zero in, non-finite out, as HLSL division).
  • quaternion.AxisAngle(float3 unitAxis, float radians): (axis sin(r/2), cos(r/2)); the axis is not normalised.
  • math.slerp(a, b, t): shortest arc (negate b when dot(a, b) < 0); normalised lerp when the quaternions are nearly parallel; t is not clamped.
  • math.swing_twist(q, unitTwistAxis, out swing, out twist): twist is q's projection onto the axis, normalised (identity when that projection vanishes, a half-turn swing); swing = mul(q, conjugate(twist)), so q == mul(swing, twist) (twist applied first).
  • math.clamp_twist(twist, unitAxis, min, max): angle 2 atan2(dot(twist.xyz, axis), twist.w) wrapped to (-pi, pi]; inside [min, max] the input is returned unchanged, otherwise AxisAngle(axis, clamp(angle, min, max)).
  • math.clamp_hinge(q, unitAxis, min, max): clamp_twist of q's twist about the axis; the swing is discarded.
  • math.clamp_cone(swing, unitAxis, maxAngle): the angle between rotate(swing, axis) and axis; inside the cone the input is returned unchanged, otherwise the rotation by maxAngle about normalize(cross(axis, rotate(swing, axis))) (at a half turn, where that cross product vanishes, about the swing's own axis).
  • Preconditions (min <= max, unit axes) are documented, not checked, as everywhere in math.

What GameCult.Animation will need from this cut (not mapped here)

  • From cut 4: mul, rotate, AxisAngle for AimChain.Pose and EffectorDirection; conjugate or inverse for frame changes; swing_twist and the clamps if its joints ever carry cone or twist limits (the presenter's per-axis boxes do not).
  • From cut 1: start-independent BoundedLeastSquares at a tight tolerance. AimChainSolver.Solve is warm-started from current; without cut 1 its answer would depend on the previous frame exactly as the allocator's did.
  • A CultMath release that carries cut 4, which org.gamecult.animation's package.json pins. That release is mapped with GameCult.Animation, so one Soul-verified release pair ships both.

Observed, not owned by these cuts

  • At the default tolerance the solver's answer on allocator-shaped problems depends on the start by up to the whole box (B6). That is documented and pinned (AnAllocatorShapedProblemIsStartIndependentAtATightTolerance); a consumer that needs one answer passes a tight tolerance, as Aetheria does.