diff --git a/crates/sketch/src/profiles.rs b/crates/sketch/src/profiles.rs index a54fc8cd..b19c35e5 100644 --- a/crates/sketch/src/profiles.rs +++ b/crates/sketch/src/profiles.rs @@ -343,7 +343,17 @@ pub fn find_profiles(sk: &Sketch) -> Vec { let on_x = x.param_on(p, tol).or_else(|| end_param(x, p, tol)); let on_y = y.param_on(p, tol).or_else(|| end_param(y, p, tol)); // Crossings at (or next to) curve ends are already cuts. - if on_x.is_none() || on_y.is_none() || ends.iter().any(|e| e.dist(p) < tol * 100.0) { + let near_end = || { + let (cx, cy) = grid.cell(p); + (-1..=1).any(|dx| { + (-1..=1).any(|dy| { + end_cells + .get(&(cx + dx, cy + dy)) + .is_some_and(|v| v.iter().any(|k| ends.get(*k).is_some_and(|e| e.dist(p) < tol * 100.0))) + }) + }) + }; + if on_x.is_none() || on_y.is_none() || near_end() { continue; } if let (Some(tx), Some(ty)) = (on_x, on_y) { @@ -611,7 +621,7 @@ pub fn find_profiles(sk: &Sketch) -> Vec { v.push(k); } } - let mut out: Vec = loops + let out: Vec = loops .iter() .enumerate() .map(|(i, (l, ids, _))| { @@ -624,8 +634,11 @@ pub fn find_profiles(sk: &Sketch) -> Vec { Profile { region, outer_curves: ids.clone(), hole_curves, area, centroid } }) .collect(); - let first_curve = |p: &Profile| p.outer_curves.iter().filter_map(|id| sk.curve_index(id)).min().unwrap_or(usize::MAX); - out.sort_by(|a, b| first_curve(a).cmp(&first_curve(b)).then(b.area.total_cmp(&a.area))); + let index: std::collections::HashMap<&str, usize> = sk.curves.iter().enumerate().map(|(i, c)| (c.id.as_str(), i)).collect(); + let first_curve = |p: &Profile| p.outer_curves.iter().filter_map(|id| index.get(id.as_str()).copied()).min().unwrap_or(usize::MAX); + let mut keyed: Vec<(usize, Profile)> = out.into_iter().map(|p| (first_curve(&p), p)).collect(); + keyed.sort_by(|a, b| a.0.cmp(&b.0).then(b.1.area.total_cmp(&a.1.area))); + let out: Vec = keyed.into_iter().map(|x| x.1).collect(); out } diff --git a/crates/sketch/src/solver.rs b/crates/sketch/src/solver.rs index 87890ef2..f9a6d2f4 100644 --- a/crates/sketch/src/solver.rs +++ b/crates/sketch/src/solver.rs @@ -444,6 +444,14 @@ enum Block<'a> { TangentAt(usize, &'a ConstraintKind, usize), } +fn copy_block<'a>(b: &Block<'a>) -> Block<'a> { + match b { + Block::Arc { c, a, b } => Block::Arc { c: *c, a: *a, b: *b }, + Block::User(i, k) => Block::User(*i, k), + Block::TangentAt(i, k, p) => Block::TangentAt(*i, k, *p), + } +} + fn block_residuals(s: &State, b: &Block, out: &mut Vec) { match b { Block::Arc { c, a, b } => { @@ -634,6 +642,53 @@ fn system(sk: &Sketch, lay: &Layout, blocks: &[(Block, Vec)], x: &mut [f6 (r, Some(j), owner) } +/// Residuals and the Jacobian restricted to `vars` (columns in that order; `col` maps a global +/// variable to its column, usize::MAX elsewhere). Blocks only touch their component's variables, +/// so a component's Jacobian is small even in a huge sketch. +fn local_system( + sk: &Sketch, + lay: &Layout, + blocks: &[(Block, Vec)], + x: &mut [f64], + col: &[usize], + ncols: usize, +) -> (Vec, Mat, Vec) { + let (r, _, owner) = system(sk, lay, blocks, x, false); + let mut j = Mat::zeros(r.len(), ncols); + let mut row = 0; + let mut plus = Vec::new(); + let mut minus = Vec::new(); + for (b, vars) in blocks { + plus.clear(); + block_residuals(&State { sk, lay, x }, b, &mut plus); + let rows = plus.len(); + for &v in vars { + let Some(c) = col.get(v).copied().filter(|c| *c < ncols) else { continue }; + let Some(x0) = x.get(v).copied() else { continue }; + let h = 1e-7 * x0.abs().max(1.0); + if let Some(xv) = x.get_mut(v) { + *xv = x0 + h; + } + plus.clear(); + block_residuals(&State { sk, lay, x }, b, &mut plus); + if let Some(xv) = x.get_mut(v) { + *xv = x0 - h; + } + minus.clear(); + block_residuals(&State { sk, lay, x }, b, &mut minus); + if let Some(xv) = x.get_mut(v) { + *xv = x0; + } + for k in 0..rows { + let d = (plus.get(k).copied().unwrap_or(0.0) - minus.get(k).copied().unwrap_or(0.0)) / (2.0 * h); + j.set(row + k, c, d); + } + } + row += rows; + } + (r, j, owner) +} + fn max_abs(v: &[f64]) -> f64 { v.iter().fold(0.0_f64, |m, x| if x.is_finite() { m.max(x.abs()) } else { f64::INFINITY }) } @@ -642,7 +697,7 @@ fn norm2(v: &[f64]) -> f64 { } /// Levenberg–Marquardt on one component (variables `vars`, blocks `blocks`). Returns iterations. -fn lm(sk: &Sketch, lay: &Layout, blocks: &[(Block, Vec)], vars: &[usize], x: &mut Vec) -> usize { +fn lm(sk: &Sketch, lay: &Layout, blocks: &[(Block, Vec)], vars: &[usize], x: &mut Vec, col: &[usize]) -> usize { let n = vars.len(); let mut lambda = 1e-3; let (mut r, _, _) = system(sk, lay, blocks, x, false); @@ -651,13 +706,13 @@ fn lm(sk: &Sketch, lay: &Layout, blocks: &[(Block, Vec)], vars: &[usize], if max_abs(&r) < TOL { return it; } - let (_, Some(j), _) = system(sk, lay, blocks, x, true) else { return it }; - // Normal equations restricted to this component's variables. + let (_, j, _) = local_system(sk, lay, blocks, x, col, n); + // Normal equations of this component. let mut a = Mat::zeros(n, n); let mut g = vec![0.0; n]; for row in 0..j.rows { let ri = r.get(row).copied().unwrap_or(0.0); - let nz: Vec<(usize, f64)> = vars.iter().enumerate().map(|(li, &v)| (li, j.get(row, v))).filter(|(_, d)| *d != 0.0).collect(); + let nz: Vec<(usize, f64)> = (0..n).map(|li| (li, j.get(row, li))).filter(|(_, d)| *d != 0.0).collect(); for &(li, di) in &nz { if let Some(gi) = g.get_mut(li) { *gi += di * ri; @@ -766,11 +821,19 @@ pub fn solve(sk: &mut Sketch) -> SolveReport { } } let mut comps: Vec<(usize, Vec, Vec)> = Vec::new(); // (root, vars, block idx) + let mut comp_of: std::collections::HashMap = std::collections::HashMap::new(); for v in 0..lay.n { let r = find(&mut parent, v); - match comps.iter_mut().find(|c| c.0 == r) { - Some(c) => c.1.push(v), - None => comps.push((r, vec![v], Vec::new())), + match comp_of.get(&r) { + Some(&k) => { + if let Some(c) = comps.get_mut(k) { + c.1.push(v); + } + } + None => { + comp_of.insert(r, comps.len()); + comps.push((r, vec![v], Vec::new())); + } } } let mut constant_blocks = Vec::new(); @@ -778,13 +841,22 @@ pub fn solve(sk: &mut Sketch) -> SolveReport { match vars.first() { Some(&v) => { let r = find(&mut parent, v); - if let Some(c) = comps.iter_mut().find(|c| c.0 == r) { + if let Some(c) = comp_of.get(&r).and_then(|k| comps.get_mut(*k)) { c.2.push(bi); } } None => constant_blocks.push(bi), } } + // Column of each variable within its component. + let mut col = vec![usize::MAX; lay.n]; + for (_, vars, _) in &comps { + for (k, v) in vars.iter().enumerate() { + if let Some(c) = col.get_mut(*v) { + *c = k; + } + } + } let mut iterations = 0; // Solve each component; `blocks` is consumed by index so rebuild per component. @@ -793,42 +865,53 @@ pub fn solve(sk: &mut Sketch) -> SolveReport { if bis.is_empty() { continue; } - let sub: Vec<(Block, Vec)> = bis - .iter() - .filter_map(|&bi| { - all_blocks.get(bi).map(|(b, v)| { - let b2 = match b { - Block::Arc { c, a, b } => Block::Arc { c: *c, a: *a, b: *b }, - Block::User(i, k) => Block::User(*i, k), - Block::TangentAt(i, k, p) => Block::TangentAt(*i, k, *p), - }; - (b2, v.clone()) - }) - }) - .collect(); - iterations += lm(sk, &lay, &sub, vars, &mut x); + let sub: Vec<(Block, Vec)> = bis.iter().filter_map(|&bi| all_blocks.get(bi).map(|(b, v)| (copy_block(b), v.clone()))).collect(); + iterations += lm(sk, &lay, &sub, vars, &mut x, &col); } - // Final residuals, Jacobian and analysis over everything. - let (r, j, owner) = system(sk, &lay, &all_blocks, &mut x, true); - let max_residual = max_abs(&r); + // Final residuals and analysis, component by component (the Jacobian is block diagonal: + // its rank is the sum of the components' ranks). + let mut max_residual = 0.0_f64; let mut failing: Vec = Vec::new(); - for (row, v) in r.iter().enumerate() { - if !(v.abs() < ACCEPT) - && let Some((Block::User(ci, _) | Block::TangentAt(ci, _, _), _)) = owner.get(row).and_then(|bi| all_blocks.get(*bi)) - && let Some(c) = sk.constraints.get(*ci) - && !failing.contains(&c.id) - { - failing.push(c.id.clone()); + let mut det = vec![false; lay.n]; + let mut rank = 0; + let mut check = |r: &[f64], owner: &[usize], blocks: &[(Block, Vec)], failing: &mut Vec| { + max_residual = max_residual.max(max_abs(r)); + for (row, v) in r.iter().enumerate() { + if !(v.abs() < ACCEPT) + && let Some((Block::User(ci, _) | Block::TangentAt(ci, _, _), _)) = owner.get(row).and_then(|bi| blocks.get(*bi)) + && let Some(c) = sk.constraints.get(*ci) + && !failing.contains(&c.id) + { + failing.push(c.id.clone()); + } + } + }; + for (_, vars, bis) in &comps { + let sub: Vec<(Block, Vec)> = bis.iter().filter_map(|&bi| all_blocks.get(bi).map(|(b, v)| (copy_block(b), v.clone()))).collect(); + if sub.is_empty() { + continue; + } + let (r, j, owner) = local_system(sk, &lay, &sub, &mut x, &col, vars.len()); + check(&r, &owner, &sub, &mut failing); + if j.rows > 0 { + let (d, rk) = determined_vars(&j, 1e-7); + rank += rk; + for (k, v) in vars.iter().enumerate() { + if let (Some(slot), Some(flag)) = (det.get_mut(*v), d.get(k)) { + *slot = *flag; + } + } } } - let _ = constant_blocks; + let consts: Vec<(Block, Vec)> = + constant_blocks.iter().filter_map(|&bi| all_blocks.get(bi).map(|(b, v)| (copy_block(b), v.clone()))).collect(); + if !consts.is_empty() { + let (r, _, owner) = system(sk, &lay, &consts, &mut x, false); + check(&r, &owner, &consts, &mut failing); + } let status = if max_residual < ACCEPT { SolveStatus::Solved } else { SolveStatus::Failed }; let use_x = if status == SolveStatus::Solved { &x } else { &x0 }; - let (det, rank) = match &j { - Some(j) if j.rows > 0 => determined_vars(j, 1e-7), - _ => (vec![false; lay.n], 0), - }; // Write back. for (i, p) in sk.points.iter_mut().enumerate() { if let Some(v) = lay.pvar.get(i).copied().flatten() { diff --git a/docs/perf.md b/docs/perf.md new file mode 100644 index 00000000..90446d77 --- /dev/null +++ b/docs/perf.md @@ -0,0 +1,36 @@ +# Performance + +## Large sketches + +Measured with `cargo test -p solvecraft-engine --release sketch_perf -- --ignored --nocapture` +(`crates/engine/src/cmd/sketch_perf_tests.rs`). The sketch holds `n` rectangles (four lines, +horizontal/vertical constraints and two length dimensions each) and `n` circles. All times are +wall-clock milliseconds on the development machine while other builds were running (load +average around 60), so treat them as upper bounds. + +| n | curves | points | constraints | solve | profiles | draw lists | snap query | add a line (command) | drag a point (command) | +|---:|---:|---:|---:|---:|---:|---:|---:|---:|---:| +| 250 | 1 250 | 1 251 | 1 500 | 15 | 30 | 12 | 0.1 | 68 | 48 | +| 500 | 2 500 | 2 501 | 3 000 | 35 | 47 | 10 | 0.1 | 123 | 100 | +| 1000 | 5 000 | 5 001 | 6 000 | 65 | 134 | 1 | 0.3 | 232 | 100 | + +"Add a line" and "drag a point" are whole commands: they copy the sketch, solve it, store it and +re-evaluate the timeline (a second solve, plus profile finding). "Draw lists" is +`view::sketch_lines`, the polylines the viewport draws. + +Before the 2026-10-08 changes the n = 1000 sketch took 1 280 ms to solve, 590 ms to find +profiles and 5.1 s per command: + +- The solver built one dense Jacobian over every variable of the sketch, for each connected + component's iterations and again for the degree-of-freedom analysis (5 000 × 10 000 entries, + then reduced row echelon form). It now builds each component's Jacobian over that component's + own variables. The sketch's Jacobian is block diagonal, so its rank is the sum of the + components' ranks, and components were already solved separately. +- Grouping variables into components searched a list per variable. It now uses a map. +- Profile finding looked at every curve end for each crossing it found, and sorted profiles with + a linear id lookup in the comparison. It now uses the grid it already builds for crossings, + and sorts on precomputed keys. + +Projection (linked reference geometry) is cached per model state, so edits to a sketch with +projected silhouettes or sections on a ~20 000-triangle model cost about 3 ms each instead of +about 80 ms.