Sketch: accept redundant but consistent constraints (#687) (#699)

A clover of same-size rectangles drawn around a shared corner was refused as over-constrained.

Constraints that remove no degree of freedom but still hold (zero residual) are now kept and shown dimmed as redundant. Conflicting ones are still refused. A redundant driving dimension is refused with "already determined by ...". Typed sizes become driven (reference) dimensions. Drawing no longer adds Horizontal/Vertical to a side that is already implied through two coincident corners. The solver reports redundant constraints using a Jacobian row-dependence test.

Co-Authored-By: Claude Opus 5.5 <[email protected]>

crates/sketch/src/lib.rs

@@ -55,8 +55,8 @@ pub use profiles::{
 };
 pub use shape::{Shape, intersections};
 pub use solver::{
-    Analysis, DrivingValues, Hold, SolveReport, SolveStatus, measure_dimension, solve, solve_after_edit, solve_after_value_edit, solve_holding,
-    solve_with,
+    Analysis, CONSISTENT, DrivingValues, Hold, SolveReport, SolveStatus, constraint_residual, determined_by, measure_dimension, solve,
+    solve_after_edit, solve_after_value_edit, solve_holding, solve_with,
 };
 pub use text::{MAX_TEXT_CHARS, text_geometry, text_geometry_in};
 pub use wire::{MAX_FIT_POINTS, MAX_WIRE_POINTS, Wire, WireSpace, fit_polyline};

crates/sketch/src/redundancy.rs

@@ -0,0 +1,171 @@
+//! Redundant constraints: ones whose equations are already implied by the others where the
+//! sketch is now (two rectangles sharing a side both say the side is vertical; a coincidence
+//! between two corners that equal sizes already put together). They remove no degree of freedom
+//! but they do not conflict either: the solver satisfies them along with the rest.
+//!
+//! The test is the Jacobian's: walking the rows in order (the curves' own relations, then the
+//! constraints in the order they were added), a row that adds no rank to the rows before it is
+//! dependent on them. A constraint none of whose rows add rank is redundant. Whether a
+//! redundant constraint is consistent is a question of its residual (see
+//! [`constraint_residual`]): zero where the other constraints hold means it is implied, not in
+//! conflict.
+
+use super::{Block, Hold, Layout, State, block_residuals, build_blocks, columns, initial_x, layout, local_system, max_abs, points_on_curves};
+use crate::linalg::{Mat, solve as lin_solve};
+use crate::model::Sketch;
+use crate::union_find::{find, union};
+
+/// A row is dependent on the rows before it when what is left of it, once their directions are
+/// taken out, is below this share of its length (the rank analysis' pivot tolerance is 1e-7 on
+/// rows of a comparable scale).
+const DEPENDENT: f64 = 1e-6;
+/// Variables in one component past which its constraints are not checked for redundancy: the
+/// check is dense, rows × rank × variables.
+pub(super) const MAX_REDUNDANCY_VARS: usize = 400;
+/// A constraint's residual (mm, or radians for angles) below which it holds.
+pub const CONSISTENT: f64 = 1e-6;
+
+/// The constraint a block stands for.
+pub(super) fn block_constraint(b: &Block) -> Option<usize> {
+    match b {
+        Block::User(ci, _) | Block::TangentAt(ci, _, _) => Some(*ci),
+        Block::Closure { .. } | Block::Arc { .. } => None,
+    }
+}
+
+fn dot(a: &[f64], b: &[f64]) -> f64 {
+    a.iter().zip(b).map(|(x, y)| x * y).sum()
+}
+
+/// Per block of `blocks` (in order, over the variables `vars`): whether every one of its rows
+/// is a combination of the rows before it. With `support`, also for the block at that index
+/// (when redundant) the earlier blocks its rows combine.
+pub(super) fn dependent_blocks(
+    sk: &Sketch,
+    lay: &Layout,
+    blocks: &[(Block, Vec<usize>)],
+    vars: &[usize],
+    x: &mut [f64],
+    support: Option<usize>,
+) -> (Vec<bool>, Vec<usize>) {
+    let mut col = vec![usize::MAX; lay.n];
+    columns(&mut col, vars, vars);
+    let (_, j, owner) = local_system(sk, lay, blocks, x, &col, vars.len());
+    let n = j.cols;
+    let row = |k: usize| -> Vec<f64> { j.data.get(k * n..(k + 1) * n).map(<[f64]>::to_vec).unwrap_or_default() };
+    // Orthonormal directions of the independent rows so far, and those rows with their block.
+    let mut q: Vec<Vec<f64>> = Vec::new();
+    let mut basis: Vec<(Vec<f64>, usize)> = Vec::new();
+    let mut has_rows = vec![false; blocks.len()];
+    let mut adds_rank = vec![false; blocks.len()];
+    let mut target_rows: Vec<Vec<f64>> = Vec::new();
+    for k in 0..j.rows {
+        let Some(&b) = owner.get(k) else { continue };
+        let v = row(k);
+        let norm = dot(&v, &v).sqrt();
+        let mut w = v.clone();
+        // Twice, for the orthogonality lost to rounding.
+        for _ in 0..2 {
+            for qi in &q {
+                let d = dot(qi, &w);
+                for (wk, qk) in w.iter_mut().zip(qi) {
+                    *wk -= d * qk;
+                }
+            }
+        }
+        let left = dot(&w, &w).sqrt();
+        if let Some(h) = has_rows.get_mut(b) {
+            *h = true;
+        }
+        if norm.is_finite() && left.is_finite() && norm > 0.0 && left > DEPENDENT * norm {
+            q.push(w.iter().map(|wk| wk / left).collect());
+            basis.push((v, b));
+            if let Some(a) = adds_rank.get_mut(b) {
+                *a = true;
+            }
+        } else if Some(b) == support {
+            target_rows.push(v);
+        }
+    }
+    let dependent: Vec<bool> = has_rows.iter().zip(&adds_rank).map(|(h, a)| *h && !*a).collect();
+    let mut by = Vec::new();
+    if let Some(t) = support
+        && dependent.get(t).copied().unwrap_or(false)
+        && !basis.is_empty()
+    {
+        // The combination of the independent rows that gives each of the target's rows: the
+        // least-squares solution of Bᵀ a = v, from (B Bᵀ) a = B v (B has full row rank).
+        let m = basis.len();
+        let mut g = Mat::zeros(m, m);
+        for (r, (br, _)) in basis.iter().enumerate() {
+            for (c, (bc, _)) in basis.iter().enumerate() {
+                g.set(r, c, dot(br, bc));
+            }
+        }
+        for v in &target_rows {
+            let rhs: Vec<f64> = basis.iter().map(|(b, _)| dot(b, v)).collect();
+            let Some(a) = lin_solve(g.clone(), rhs) else { continue };
+            let big = a.iter().fold(0.0_f64, |m, x| m.max(x.abs()));
+            for (ak, (_, b)) in a.iter().zip(&basis) {
+                if ak.abs() > 1e-6 * big && !by.contains(b) {
+                    by.push(*b);
+                }
+            }
+        }
+        by.sort_unstable();
+    }
+    (dependent, by)
+}
+
+/// How far the constraint at index `ci` is from holding where the sketch is now: its largest
+/// residual (mm, or radians for angles). Zero for a constraint that holds, so a constraint that
+/// removes no degree of freedom but holds is redundant, not conflicting. `None` for a driven
+/// dimension (it constrains nothing) or an unknown index.
+pub fn constraint_residual(sk: &Sketch, ci: usize) -> Option<f64> {
+    let c = sk.constraints.get(ci).filter(|c| !c.driven)?;
+    let fixed = Layout { pvar: vec![None; sk.points.len()], rvar: vec![None; sk.curves.len()], n: 0, points: 0 };
+    let on_curve = points_on_curves(sk);
+    let b = match super::touch_point(sk, &c.kind, &on_curve) {
+        Some(p) => Block::TangentAt(ci, &c.kind, p),
+        None => Block::User(ci, &c.kind),
+    };
+    let mut r = Vec::new();
+    block_residuals(&State { sk, lay: &fixed, x: &[] }, &b, &mut r);
+    Some(max_abs(&r))
+}
+
+/// The constraints (indices) that already determine what the constraint at index `ci` asks
+/// for, where the sketch is now: the ones whose equations its own are a combination of. Empty
+/// when it is not redundant (it removes a degree of freedom), or when its component is too
+/// large to check.
+pub fn determined_by(sk: &Sketch, ci: usize) -> Vec<usize> {
+    let lay = layout(sk, &Hold::default());
+    let mut x = initial_x(sk, &lay);
+    let on_curve = points_on_curves(sk);
+    let (blocks, _) = build_blocks(sk, &lay, &on_curve);
+    let Some(t) = blocks.iter().position(|(b, _)| block_constraint(b) == Some(ci)) else { return Vec::new() };
+    let mut parent: Vec<usize> = (0..lay.n).collect();
+    for (_, vars) in &blocks {
+        if let Some(&first) = vars.first() {
+            for &v in vars.iter().skip(1) {
+                union(&mut parent, first, v);
+            }
+        }
+    }
+    let Some(root) = blocks.get(t).and_then(|(_, v)| v.first().copied()).map(|v| find(&mut parent, v)) else { return Vec::new() };
+    let vars: Vec<usize> = (0..lay.n).filter(|v| find(&mut parent, *v) == root).collect();
+    if vars.len() > MAX_REDUNDANCY_VARS {
+        return Vec::new();
+    }
+    // The component's other blocks in order, then the constraint asked about.
+    let mut order: Vec<usize> = (0..blocks.len())
+        .filter(|bi| *bi != t && blocks.get(*bi).and_then(|(_, v)| v.first()).is_some_and(|v| find(&mut parent, *v) == root))
+        .collect();
+    order.push(t);
+    let sub: Vec<(Block, Vec<usize>)> = order.iter().filter_map(|bi| blocks.get(*bi).map(|(b, v)| (super::copy_block(b), v.clone()))).collect();
+    let (_, by) = dependent_blocks(sk, &lay, &sub, &vars, &mut x, Some(sub.len().saturating_sub(1)));
+    let mut out: Vec<usize> = by.iter().filter_map(|k| sub.get(*k).and_then(|(b, _)| block_constraint(b))).collect();
+    out.sort_unstable();
+    out.dedup();
+    out
+}

crates/sketch/src/solver.rs

@@ -21,6 +21,10 @@ use crate::linalg::{Mat, determined_vars, solve as lin_solve};
 use crate::model::{ConstraintKind, CurveKind, Sketch};
 use crate::union_find::{find, union};
 
+#[path = "redundancy.rs"]
+mod redundancy;
+pub use redundancy::{CONSISTENT, constraint_residual, determined_by};
+
 const MAX_ITERS: usize = 500;
 /// Variables in one connected component. Each component is solved as a dense system: an n by n
 /// matrix for the Levenberg-Marquardt step, rows by n for the Jacobian. 4,000 variables is a
@@ -65,6 +69,12 @@ pub struct SolveReport {
     pub point_determined: Vec<bool>,
     /// Per curve: all of its defining values determined.
     pub curve_determined: Vec<bool>,
+    /// Constraints that remove no degree of freedom: everything they ask for is already implied
+    /// by the constraints before them (a second "vertical" on a side two rectangles share).
+    /// They hold, so they are not a conflict; they are shown as redundant. Only from a full
+    /// analysis, and only for components of at most a few hundred variables.
+    #[serde(default)]
+    pub redundant: Vec<String>,
 }
 
 impl SolveReport {
@@ -1298,6 +1308,60 @@ impl Sketch {
     }
 }
 
+/// The variables' values where the sketch is now.
+fn initial_x(sk: &Sketch, lay: &Layout) -> Vec<f64> {
+    let mut x = vec![0.0; lay.n];
+    for (i, p) in sk.points.iter().enumerate() {
+        if let Some(v) = lay.pvar.get(i).copied().flatten() {
+            set_point(&mut x, v, p.pos);
+        }
+    }
+    for (i, c) in sk.curves.iter().enumerate() {
+        if let (CurveKind::Circle { r, .. } | CurveKind::Ellipse { r, .. }, Some(v)) = (&c.kind, lay.rvar.get(i).copied().flatten())
+            && let Some(s) = x.get_mut(v)
+        {
+            *s = *r;
+        }
+    }
+    x
+}
+
+/// The solve's blocks in order (the curves' own relations, then the driving constraints in
+/// constraint order) with the variables each reads, and per block whether it reads a constant
+/// point.
+fn build_blocks<'a>(sk: &'a Sketch, lay: &Layout, on_curve: &[(usize, usize)]) -> (Vec<(Block<'a>, Vec<usize>)>, Vec<bool>) {
+    let mut raw = Vec::new();
+    for c in &sk.curves {
+        if let CurveKind::Spline { pts, periodic: true, .. } = &c.kind
+            && let (Some(a), Some(b)) = (pts.first(), pts.last())
+            && a != b
+        {
+            raw.push(Block::Closure { a: *a, b: *b });
+        }
+        if let CurveKind::Arc { c, a, b } = c.kind {
+            raw.push(Block::Arc { c, a, b });
+        }
+    }
+    for (i, c) in sk.constraints.iter().enumerate() {
+        if c.driven {
+            continue;
+        }
+        raw.push(match touch_point(sk, &c.kind, on_curve) {
+            Some(p) => Block::TangentAt(i, &c.kind, p),
+            None => Block::User(i, &c.kind),
+        });
+    }
+    let mut blocks: Vec<(Block, Vec<usize>)> = Vec::with_capacity(raw.len());
+    let mut reads_constant = Vec::with_capacity(raw.len());
+    for bl in raw {
+        let (vars, constant) = block_vars(sk, lay, &bl);
+        blocks.push((bl, vars));
+        reads_constant.push(constant);
+    }
+
+    (blocks, reads_constant)
+}
+
 #[cfg(not(test))]
 fn solve_from(
     sk: &mut Sketch,
@@ -1721,19 +1785,7 @@ fn solve_inner(
         }
     }
     let lay = layout(sk, hold);
-    let mut x = vec![0.0; lay.n];
-    for (i, p) in sk.points.iter().enumerate() {
-        if let Some(v) = lay.pvar.get(i).copied().flatten() {
-            set_point(&mut x, v, p.pos);
-        }
-    }
-    for (i, c) in sk.curves.iter().enumerate() {
-        if let (CurveKind::Circle { r, .. } | CurveKind::Ellipse { r, .. }, Some(v)) = (&c.kind, lay.rvar.get(i).copied().flatten())
-            && let Some(s) = x.get_mut(v)
-        {
-            *s = *r;
-        }
-    }
+    let mut x = initial_x(sk, &lay);
     let x0 = x.clone();
     let on_curve = points_on_curves(sk);
     // Start from the constructed positions of the points the constraints pin down (see
@@ -1785,34 +1837,7 @@ fn solve_inner(
     nudge_line_angles(sk, &lay, &mut x);
 
     // The blocks, and per block whether it reads a constant point.
-    let mut raw = Vec::new();
-    for c in &sk.curves {
-        if let CurveKind::Spline { pts, periodic: true, .. } = &c.kind
-            && let (Some(a), Some(b)) = (pts.first(), pts.last())
-            && a != b
-        {
-            raw.push(Block::Closure { a: *a, b: *b });
-        }
-        if let CurveKind::Arc { c, a, b } = c.kind {
-            raw.push(Block::Arc { c, a, b });
-        }
-    }
-    for (i, c) in sk.constraints.iter().enumerate() {
-        if c.driven {
-            continue;
-        }
-        raw.push(match touch_point(sk, &c.kind, &on_curve) {
-            Some(p) => Block::TangentAt(i, &c.kind, p),
-            None => Block::User(i, &c.kind),
-        });
-    }
-    let mut blocks: Vec<(Block, Vec<usize>)> = Vec::with_capacity(raw.len());
-    let mut reads_constant = Vec::with_capacity(raw.len());
-    for bl in raw {
-        let (vars, constant) = block_vars(sk, &lay, &bl);
-        blocks.push((bl, vars));
-        reads_constant.push(constant);
-    }
+    let (blocks, reads_constant) = build_blocks(sk, &lay, &on_curve);
 
     carry_with_seeded(&mut x, &x0, &constructed, &blocks, &lay);
     // The plain start, for the variables whose start the seed (or carrying along it) moved.
@@ -1981,6 +2006,7 @@ fn solve_inner(
     let mut failing: Vec<String> = Vec::new();
     let mut det = vec![false; lay.n];
     let mut rank = 0;
+    let mut redundant: Vec<String> = Vec::new();
     // A component too large to solve is a failure, not an allocation.
     if comps.iter().any(|(_, vars, bis)| !bis.is_empty() && oversized(vars)) {
         max_residual = f64::MAX;
@@ -2023,6 +2049,17 @@ fn solve_inner(
         let (rest_det, rest_rank) = if j.rows > 0 { determined_vars(&j, 1e-7) } else { (vec![false; rest.len()], 0) };
         let comp_rank = pin_flags.iter().filter(|p| **p).count() + rest_rank;
         rank += comp_rank;
+        // More rows than rank: some constraint may be implied by the others.
+        if r.len() > comp_rank && vars.len() <= redundancy::MAX_REDUNDANCY_VARS {
+            let (dependent, _) = redundancy::dependent_blocks(sk, &lay, &sub, vars, &mut x, None);
+            for ((b, _), d) in sub.iter().zip(dependent) {
+                if d && let Some(c) = redundancy::block_constraint(b).and_then(|ci| sk.constraints.get(ci))
+                    && !redundant.contains(&c.id)
+                {
+                    redundant.push(c.id.clone());
+                }
+            }
+        }
         for (&v, pinned) in vars.iter().zip(&pin_flags) {
             let flag = match rest.binary_search(&v) {
                 Ok(k) => rest_det.get(k).copied().unwrap_or(false),
@@ -2171,6 +2208,7 @@ fn solve_inner(
         continuations,
         point_determined,
         curve_determined,
+        redundant,
     }
 }