create/docs/src/solver/diagnostics.md

Diagnostics

The solver provides three levels of constraint analysis: system-wide DOF counting, per-entity DOF decomposition, and overconstrained/conflicting constraint detection.

DOF counting

Source: mods/solver/kindred_solver/dof.py

Degrees of freedom are computed from the Jacobian rank:

DOF = n_free_params - rank(J)

Where n_free_params is the number of non-fixed parameters and rank(J) is the numerical rank of the Jacobian evaluated at current parameter values (SVD with tolerance 1e-8).

A well-constrained assembly has DOF = 0 (exactly enough constraints to determine all positions). Positive DOF means underconstrained (parts can still move). Negative DOF is not possible with this formulation -- instead, rank deficiency in an overdetermined system indicates redundant constraints.

The DOF value is reported in SolveResult.dof after every solve.

Per-entity DOF

Source: mods/solver/kindred_solver/diagnostics.py

per_entity_dof() breaks down the DOF count per body, identifying which motions remain free for each part:

1. Build the full Jacobian
2. For each non-grounded body, extract the 7 columns corresponding to its parameters
3. Compute SVD of the sub-matrix; rank = number of constrained directions
4. remaining_dof = 7 - rank (includes the quaternion normalization constraint counted in the rank)
5. Classify null-space vectors as free motions by analyzing their translation vs. rotation components:
   • Pure translation: >80% of the null vector's energy is in tx, ty, tz components
   • Pure rotation: >80% of the energy is in qw, qx, qy, qz components
   • Helical: mixed

Returns a list of EntityDOF dataclasses:

@dataclass
class EntityDOF:
    entity_id: str
    remaining_dof: int
    free_motions: list[str]  # e.g., ["rotation about Z", "translation along X"]

Overconstrained detection

Source: mods/solver/kindred_solver/diagnostics.py

find_overconstrained() identifies redundant and conflicting constraints when the system is overconstrained (Jacobian is rank-deficient). It runs automatically when solve() fails to converge.

Algorithm

Following the approach used by SolvSpace:

1. Check rank. Build the full Jacobian J, compute its rank via SVD. If rank == n_residuals, the system is not overconstrained -- return empty.

2. Find redundant constraints. For each constraint, temporarily remove its rows from J and re-check rank. If the rank is preserved, the constraint is redundant (removing it doesn't change the system's effective equations).

3. Distinguish conflicting from merely redundant. Compute the left null space of J (columns of U beyond the rank). Project the residual vector onto this null space:

   null_residual = U_null^T @ r
   residual_projection = U_null @ null_residual
   

   If a redundant constraint's residuals have significant projection onto the null space, it is conflicting -- it's both redundant and unsatisfied, meaning it contradicts other constraints.

Diagnostic output

Returns ConstraintDiag dataclasses:

@dataclass
class ConstraintDiag:
    constraint_index: int
    kind: str           # "redundant" or "conflicting"
    detail: str         # Human-readable explanation

These are converted to kcsolve.ConstraintDiagnostic objects in the IKCSolver bridge:

ConstraintDiag.kind	kcsolve.DiagnosticKind
"redundant"	Redundant
"conflicting"	Conflicting

Example

Two Fixed joints between the same pair of parts:

• Joint A: 6 residuals (3 position + 3 orientation)
• Joint B: 6 residuals (same as Joint A)

Jacobian rank = 6 (Joint B's rows are linearly dependent on Joint A's). Both joints are detected as redundant. If the joints specify different relative positions, both are also flagged as conflicting.

Solution preferences

Source: mods/solver/kindred_solver/preference.py

Solution preferences guide the solver toward physically intuitive solutions when multiple valid configurations exist.

Minimum-movement weighting

The weight vector scales the Newton step to prefer solutions near the initial configuration. Translation parameters get weight 1.0, quaternion parameters get weight (180/pi)^2 ~ 3283. This makes a 1-radian rotation equally "expensive" as a ~57-unit translation.

The weighted minimum-norm step is:

J_scaled = J @ diag(W^{-1/2})
dx_scaled = lstsq(J_scaled, -r)
dx = dx_scaled * W^{-1/2}

This produces the minimum-norm solution in the weighted parameter space, biasing toward small movements.

Half-space tracking

Described in detail in Solving Algorithms: Half-space tracking. Preserves the initial configuration's "branch" for constraints with multiple valid solutions by detecting and correcting branch crossings during iteration.