arael_sketch_solver/
lib.rs

1//! 2D parametric constraint-based sketch solver.
2//!
3//! Work in progress -- a parametric CAD sketching tool built on the arael
4//! optimization framework. Draw geometry, apply constraints, and the
5//! solver keeps everything consistent in real time.
6//!
7//! # Entities
8//!
9//! Three entity types, each owning its own parameters:
10//!
11//! - [`Point`] -- 2D position (x, y)
12//! - [`Line`] -- two endpoints p1, p2 (4 params), plus optional
13//!   horizontal/vertical/length constraints
14//! - [`Arc`] -- center, radius, start/end angle (5 params for arcs,
15//!   3 for circles where angles are fixed)
16//!
17//! Shared geometry (e.g. two lines meeting at a point) is enforced via
18//! coincident constraints, not shared references.
19//!
20//! # Constraints
21//!
22//! Over 40 constraint types including: coincident (point-point, line-line,
23//! point-on-line, point-on-arc), parallel, perpendicular, tangent,
24//! equal length/radius, distance, horizontal/vertical distance, and more.
25//! All constraints are symbolically differentiated at compile time.
26//!
27//! # Solving
28//!
29//! The [`Sketch::solve()`] method runs Levenberg-Marquardt optimization with
30//! drift regularization. Locking: use `Param::fixed(value)` to pin
31//! individual parameters. The solver skips fixed params entirely.
32
33pub mod dimensions;
34pub use dimensions::*;
35pub mod symbol_bag;
36pub use symbol_bag::SymbolBag;
37pub mod expr_constraint;
38pub use expr_constraint::ExpressionConstraint;
39pub mod blocker;
40pub use blocker::{BlockerReport, analyze as analyze_blockers};
41
42use arael::model::{CrossBlock, JacobianModel, Model, Param, SelfBlock, TripletBlock};
43
44const TIMING_DEBUG: bool = false;
45use arael::vect::vect2d;
46use arael::refs::{Ref, Arena};
47
48// Entity and constraint types must share one module scope because the
49// #[arael::model] macro emits private `_PARAM_COUNT` constants that
50// CrossBlock<A, B> expansions need to reference.
51include!("entities.rs");
52include!("constraints.rs");
53
54// ---------------------------------------------------------------------------
55// DOF analysis
56// ---------------------------------------------------------------------------
57
58/// Snapshot returned by `Sketch::add_drag_auto_anchors` and consumed
59/// by `Sketch::remove_drag_auto_anchors` to roll auto-anchors back at
60/// drag-end. Cloneable so the GUI can also roll anchors back from a
61/// serialized clone of the live sketch.
62#[derive(Clone)]
63pub struct DragAutoAnchorState {
64    helper_points: std::vec::Vec<arael::refs::Ref<Point>>,
65    coincident_lp1_len: usize,
66    coincident_lp2_len: usize,
67    coincident_arc_start_len: usize,
68    coincident_arc_end_len: usize,
69}
70
71/// Result of DOF (degrees of freedom) analysis.
72pub struct DofResult {
73    /// Number of unconstrained degrees of freedom.
74    pub dof: usize,
75    /// Parameter names indexed by param index (only filled when analyze=true).
76    pub param_names: Vec<String>,
77    /// Eigenvalues from Hessian decomposition (only filled when analyze=true).
78    pub eigenvalues: Vec<f64>,
79    /// Eigenvectors, one per eigenvalue (only filled when analyze=true).
80    pub eigenvectors: Vec<Vec<f64>>,
81}
82
83// ---------------------------------------------------------------------------
84// Root
85// ---------------------------------------------------------------------------
86
87#[derive(serde::Serialize, serde::Deserialize)]
88#[arael::model]
89#[arael(root, extended, jacobian)]
90pub struct Sketch {
91    pub points: Arena<Point>,
92    pub lines: Arena<Line>,
93    pub arcs: Arena<Arc>,
94    // Solver parameters
95    pub drift_isigma: f64,
96    pub constraint_isigma: f64,
97    /// Minimum length threshold for soft Heaviside penalties on line length,
98    /// arc radius, and tangent projection.
99    #[serde(default = "default_min_length")]
100    pub min_length: f64,
101    #[serde(default)]
102    pub verbose: bool,
103    // Auto-naming counters
104    pub next_point_id: u32,
105    pub next_line_id: u32,
106    pub next_arc_id: u32,
107    // Cross-constraint collections
108    pub coincident_pp: std::vec::Vec<CoincidentPP>,
109    pub coincident_lp1: std::vec::Vec<CoincidentLP1>,
110    pub coincident_lp2: std::vec::Vec<CoincidentLP2>,
111    pub coincident_ll11: std::vec::Vec<CoincidentLL11>,
112    pub coincident_ll12: std::vec::Vec<CoincidentLL12>,
113    pub coincident_ll21: std::vec::Vec<CoincidentLL21>,
114    pub coincident_ll22: std::vec::Vec<CoincidentLL22>,
115    pub distance_pp: std::vec::Vec<DistancePP>,
116    pub hdistance_pp: std::vec::Vec<HorizontalDistancePP>,
117    pub vdistance_pp: std::vec::Vec<VerticalDistancePP>,
118    pub point_on_line: std::vec::Vec<PointOnLine>,
119    pub midpoint: std::vec::Vec<MidpointConstraint>,
120    pub midpoint_lp1: std::vec::Vec<MidpointLP1>,
121    pub midpoint_lp2: std::vec::Vec<MidpointLP2>,
122    pub midpoint_arc_start: std::vec::Vec<MidpointArcStart>,
123    pub midpoint_arc_end: std::vec::Vec<MidpointArcEnd>,
124    #[serde(default)]
125    pub midpoint_arc_point: std::vec::Vec<MidpointArcPoint>,
126    #[serde(default)]
127    pub midpoint_lp1_arc: std::vec::Vec<MidpointLP1Arc>,
128    #[serde(default)]
129    pub midpoint_lp2_arc: std::vec::Vec<MidpointLP2Arc>,
130    #[serde(default)]
131    pub midpoint_arc_start_arc: std::vec::Vec<MidpointArcStartArc>,
132    #[serde(default)]
133    pub midpoint_arc_end_arc: std::vec::Vec<MidpointArcEndArc>,
134    pub point_on_arc: std::vec::Vec<PointOnArc>,
135    pub parallel: std::vec::Vec<Parallel>,
136    pub perpendicular: std::vec::Vec<Perpendicular>,
137    #[serde(default)]
138    pub arc_line_parallel: std::vec::Vec<ArcLineParallel>,
139    #[serde(default)]
140    pub arc_arc_parallel: std::vec::Vec<ArcArcParallel>,
141    pub collinear: std::vec::Vec<Collinear>,
142    pub equal_length: std::vec::Vec<EqualLength>,
143    pub angle: std::vec::Vec<AngleConstraint>,
144    pub tangent_la: std::vec::Vec<TangentLA>,
145    pub concentric: std::vec::Vec<Concentric>,
146    pub equal_radius: std::vec::Vec<EqualRadius>,
147    pub tangent_aa: std::vec::Vec<TangentAA>,
148    pub symmetry_ll: std::vec::Vec<SymmetryLL>,
149    #[serde(default)]
150    pub symmetry_pp: std::vec::Vec<SymmetryPP>,
151    #[serde(default)]
152    pub symmetry_aa: std::vec::Vec<SymmetryAA>,
153    pub distance_pl: std::vec::Vec<DistancePL>,
154    pub distance_lp1l: std::vec::Vec<DistanceLP1L>,
155    pub distance_lp2l: std::vec::Vec<DistanceLP2L>,
156    pub distance_arc_center_l: std::vec::Vec<DistanceArcCenterL>,
157    pub distance_arc_start_l: std::vec::Vec<DistanceArcStartL>,
158    pub distance_arc_end_l: std::vec::Vec<DistanceArcEndL>,
159    pub line_p1_on_line: std::vec::Vec<LineP1OnLine>,
160    pub line_p2_on_line: std::vec::Vec<LineP2OnLine>,
161    pub coincident_arc_center: std::vec::Vec<CoincidentArcCenter>,
162    pub coincident_arc_start: std::vec::Vec<CoincidentArcStart>,
163    pub coincident_arc_end: std::vec::Vec<CoincidentArcEnd>,
164    // Line endpoint <-> Arc point
165    pub coincident_lp1_arc_center: std::vec::Vec<CoincidentLP1ArcCenter>,
166    pub coincident_lp2_arc_center: std::vec::Vec<CoincidentLP2ArcCenter>,
167    pub coincident_lp1_arc_start: std::vec::Vec<CoincidentLP1ArcStart>,
168    pub coincident_lp2_arc_start: std::vec::Vec<CoincidentLP2ArcStart>,
169    pub coincident_lp1_arc_end: std::vec::Vec<CoincidentLP1ArcEnd>,
170    pub coincident_lp2_arc_end: std::vec::Vec<CoincidentLP2ArcEnd>,
171    // Arc-Arc endpoint
172    pub coincident_arc_center_start: std::vec::Vec<CoincidentArcCenterStart>,
173    pub coincident_arc_center_end: std::vec::Vec<CoincidentArcCenterEnd>,
174    pub coincident_arc_start_center: std::vec::Vec<CoincidentArcStartCenter>,
175    pub coincident_arc_end_center: std::vec::Vec<CoincidentArcEndCenter>,
176    pub coincident_arc_start_start: std::vec::Vec<CoincidentArcStartStart>,
177    pub coincident_arc_start_end: std::vec::Vec<CoincidentArcStartEnd>,
178    pub coincident_arc_end_start: std::vec::Vec<CoincidentArcEndStart>,
179    pub coincident_arc_end_end: std::vec::Vec<CoincidentArcEndEnd>,
180    pub line_p1_on_arc: std::vec::Vec<LineP1OnArc>,
181    pub line_p2_on_arc: std::vec::Vec<LineP2OnArc>,
182    pub distance_ll11: std::vec::Vec<DistanceLL11>,
183    pub distance_ll12: std::vec::Vec<DistanceLL12>,
184    pub distance_ll21: std::vec::Vec<DistanceLL21>,
185    pub distance_ll22: std::vec::Vec<DistanceLL22>,
186    pub distance_lp1: std::vec::Vec<DistanceLP1>,
187    pub distance_lp2: std::vec::Vec<DistanceLP2>,
188    pub distance_arc_center_p: std::vec::Vec<DistanceArcCenterP>,
189    pub distance_arc_start_p: std::vec::Vec<DistanceArcStartP>,
190    pub distance_arc_end_p: std::vec::Vec<DistanceArcEndP>,
191    pub distance_arc_center_l1: std::vec::Vec<DistanceArcCenterL1>,
192    pub distance_arc_center_l2: std::vec::Vec<DistanceArcCenterL2>,
193    pub distance_arc_start_l1: std::vec::Vec<DistanceArcStartL1>,
194    pub distance_arc_start_l2: std::vec::Vec<DistanceArcStartL2>,
195    pub distance_arc_end_l1: std::vec::Vec<DistanceArcEndL1>,
196    pub distance_arc_end_l2: std::vec::Vec<DistanceArcEndL2>,
197    pub distance_aa_ce_ce: std::vec::Vec<DistanceAACeCe>,
198    pub distance_aa_ce_s: std::vec::Vec<DistanceAACeS>,
199    pub distance_aa_ce_e: std::vec::Vec<DistanceAACeE>,
200    pub distance_aa_s_ce: std::vec::Vec<DistanceAASCe>,
201    pub distance_aa_s_s: std::vec::Vec<DistanceAASS>,
202    pub distance_aa_s_e: std::vec::Vec<DistanceAASE>,
203    pub distance_aa_e_ce: std::vec::Vec<DistanceAAECe>,
204    pub distance_aa_e_s: std::vec::Vec<DistanceAAES>,
205    pub distance_aa_e_e: std::vec::Vec<DistanceAAEE>,
206    /// Radial distance between two concentric arcs/circles.
207    #[serde(default)] pub distance_concentric: std::vec::Vec<DistanceConcentric>,
208    // Axis distance (horizontal/vertical unified)
209    #[serde(default)] pub axis_distance_ll11: std::vec::Vec<AxisDistanceLL11>,
210    #[serde(default)] pub axis_distance_ll12: std::vec::Vec<AxisDistanceLL12>,
211    #[serde(default)] pub axis_distance_ll21: std::vec::Vec<AxisDistanceLL21>,
212    #[serde(default)] pub axis_distance_ll22: std::vec::Vec<AxisDistanceLL22>,
213    #[serde(default)] pub axis_distance_lp1: std::vec::Vec<AxisDistanceLP1>,
214    #[serde(default)] pub axis_distance_lp2: std::vec::Vec<AxisDistanceLP2>,
215    #[serde(default)] pub axis_distance_arc_center_p: std::vec::Vec<AxisDistanceArcCenterP>,
216    #[serde(default)] pub axis_distance_arc_start_p: std::vec::Vec<AxisDistanceArcStartP>,
217    #[serde(default)] pub axis_distance_arc_end_p: std::vec::Vec<AxisDistanceArcEndP>,
218    #[serde(default)] pub axis_distance_arc_center_l1: std::vec::Vec<AxisDistanceArcCenterL1>,
219    #[serde(default)] pub axis_distance_arc_center_l2: std::vec::Vec<AxisDistanceArcCenterL2>,
220    #[serde(default)] pub axis_distance_arc_start_l1: std::vec::Vec<AxisDistanceArcStartL1>,
221    #[serde(default)] pub axis_distance_arc_start_l2: std::vec::Vec<AxisDistanceArcStartL2>,
222    #[serde(default)] pub axis_distance_arc_end_l1: std::vec::Vec<AxisDistanceArcEndL1>,
223    #[serde(default)] pub axis_distance_arc_end_l2: std::vec::Vec<AxisDistanceArcEndL2>,
224    #[serde(default)] pub axis_distance_aa_ce_ce: std::vec::Vec<AxisDistanceAACeCe>,
225    #[serde(default)] pub axis_distance_aa_ce_s: std::vec::Vec<AxisDistanceAACeS>,
226    #[serde(default)] pub axis_distance_aa_ce_e: std::vec::Vec<AxisDistanceAACeE>,
227    #[serde(default)] pub axis_distance_aa_s_ce: std::vec::Vec<AxisDistanceAASCe>,
228    #[serde(default)] pub axis_distance_aa_s_s: std::vec::Vec<AxisDistanceAASS>,
229    #[serde(default)] pub axis_distance_aa_s_e: std::vec::Vec<AxisDistanceAASE>,
230    #[serde(default)] pub axis_distance_aa_e_ce: std::vec::Vec<AxisDistanceAAECe>,
231    #[serde(default)] pub axis_distance_aa_e_s: std::vec::Vec<AxisDistanceAAES>,
232    #[serde(default)] pub axis_distance_aa_e_e: std::vec::Vec<AxisDistanceAAEE>,
233    // Dimension annotations
234    #[arael(skip)]
235    pub dimensions: std::vec::Vec<Dimension>,
236    pub next_dimension_id: u32,
237    // Next auto-assigned numeric constraint id (C<nid>). 0 is reserved
238    // as the "unassigned" sentinel picked up by assign_constraint_names().
239    #[serde(default = "default_next_constraint_id")]
240    pub next_constraint_id: u32,
241    // User-defined parameters
242    #[arael(skip)]
243    #[serde(default)]
244    pub user_params: std::vec::Vec<UserParam>,
245    // Expression constraints (parametric dimensions)
246    #[arael(skip)]
247    #[serde(skip)]
248    pub expr_constraints: std::vec::Vec<ExpressionConstraint>,
249    #[arael(skip)]
250    #[serde(skip)]
251    symbol_bag: Option<SymbolBag>,
252    // Shared TripletBlock for all expression constraints
253    #[serde(skip)]
254    pub expr_hb: TripletBlock<f64>,
255    // Cached DOF count — set by compute_dof(), cleared on structural mutation
256    #[arael(skip)]
257    #[serde(skip)]
258    pub cached_dof: Option<usize>,
259}
260
261fn default_min_length() -> f64 { 0.0001 }
262fn default_next_constraint_id() -> u32 { 1 }
263
264/// Format a synthetic constraint name for a flag-style constraint on a
265/// named entity. Pattern: "C<entity><flag>", e.g. "CL0H", "CL3V".
266pub fn format_flag_name(entity_name: &str, flag: char) -> String {
267    format!("C{}{}", entity_name, flag)
268}
269
270/// Parse a synthetic constraint name into its entity-name + flag-char
271/// components. Returns None for names that don't match "C<entity>F"
272/// where the first char is 'C' and the final char is an uppercase
273/// ASCII flag tag ('H' or 'V' today).
274pub fn parse_flag_name(token: &str) -> Option<(String, char)> {
275    let bytes = token.as_bytes();
276    if bytes.len() < 3 || bytes[0] != b'C' { return None; }
277    let flag = bytes[bytes.len() - 1] as char;
278    if !matches!(flag, 'H' | 'V') { return None; }
279    let entity = &token[1..token.len() - 1];
280    if entity.is_empty() { return None; }
281    Some((entity.to_string(), flag))
282}
283
284// ---------------------------------------------------------------------------
285// Helpers
286// ---------------------------------------------------------------------------
287
288/// Locate the cid of the axis-distance constraint (horizontal=true
289/// for HDistance, false for VDistance) that matches the endpoint
290/// pair. Mirrors the add path in
291/// `arael-sketch/src/actions.rs::push_axis_distance`.
292fn find_axis_cid(
293    s: &Sketch,
294    a: &crate::dimensions::DimensionEndpoint,
295    b: &crate::dimensions::DimensionEndpoint,
296    horizontal: bool,
297) -> Option<u32> {
298    use crate::dimensions::DimensionEndpoint::*;
299    match (a, b) {
300        (Point(pa), Point(pb)) => {
301            if horizontal {
302                s.hdistance_pp.iter().find(|c| c.a == *pa && c.b == *pb).map(|c| c.cid)
303            } else {
304                s.vdistance_pp.iter().find(|c| c.a == *pa && c.b == *pb).map(|c| c.cid)
305            }
306        }
307        (LineP1(la), LineP1(lb)) => s.axis_distance_ll11.iter().find(|c| c.a == *la && c.b == *lb && c.horizontal == horizontal).map(|c| c.cid),
308        (LineP1(la), LineP2(lb)) => s.axis_distance_ll12.iter().find(|c| c.a == *la && c.b == *lb && c.horizontal == horizontal).map(|c| c.cid),
309        (LineP2(la), LineP1(lb)) => s.axis_distance_ll21.iter().find(|c| c.a == *la && c.b == *lb && c.horizontal == horizontal).map(|c| c.cid),
310        (LineP2(la), LineP2(lb)) => s.axis_distance_ll22.iter().find(|c| c.a == *la && c.b == *lb && c.horizontal == horizontal).map(|c| c.cid),
311        (LineP1(l), Point(p)) | (Point(p), LineP1(l)) =>
312            s.axis_distance_lp1.iter().find(|c| c.line == *l && c.point == *p && c.horizontal == horizontal).map(|c| c.cid),
313        (LineP2(l), Point(p)) | (Point(p), LineP2(l)) =>
314            s.axis_distance_lp2.iter().find(|c| c.line == *l && c.point == *p && c.horizontal == horizontal).map(|c| c.cid),
315        (ArcCenter(ar), Point(p)) | (Point(p), ArcCenter(ar)) =>
316            s.axis_distance_arc_center_p.iter().find(|c| c.arc == *ar && c.point == *p && c.horizontal == horizontal).map(|c| c.cid),
317        (ArcStart(ar), Point(p)) | (Point(p), ArcStart(ar)) =>
318            s.axis_distance_arc_start_p.iter().find(|c| c.arc == *ar && c.point == *p && c.horizontal == horizontal).map(|c| c.cid),
319        (ArcEnd(ar), Point(p)) | (Point(p), ArcEnd(ar)) =>
320            s.axis_distance_arc_end_p.iter().find(|c| c.arc == *ar && c.point == *p && c.horizontal == horizontal).map(|c| c.cid),
321        (ArcCenter(ar), LineP1(l)) | (LineP1(l), ArcCenter(ar)) =>
322            s.axis_distance_arc_center_l1.iter().find(|c| c.arc == *ar && c.line == *l && c.horizontal == horizontal).map(|c| c.cid),
323        (ArcCenter(ar), LineP2(l)) | (LineP2(l), ArcCenter(ar)) =>
324            s.axis_distance_arc_center_l2.iter().find(|c| c.arc == *ar && c.line == *l && c.horizontal == horizontal).map(|c| c.cid),
325        (ArcStart(ar), LineP1(l)) | (LineP1(l), ArcStart(ar)) =>
326            s.axis_distance_arc_start_l1.iter().find(|c| c.arc == *ar && c.line == *l && c.horizontal == horizontal).map(|c| c.cid),
327        (ArcStart(ar), LineP2(l)) | (LineP2(l), ArcStart(ar)) =>
328            s.axis_distance_arc_start_l2.iter().find(|c| c.arc == *ar && c.line == *l && c.horizontal == horizontal).map(|c| c.cid),
329        (ArcEnd(ar), LineP1(l)) | (LineP1(l), ArcEnd(ar)) =>
330            s.axis_distance_arc_end_l1.iter().find(|c| c.arc == *ar && c.line == *l && c.horizontal == horizontal).map(|c| c.cid),
331        (ArcEnd(ar), LineP2(l)) | (LineP2(l), ArcEnd(ar)) =>
332            s.axis_distance_arc_end_l2.iter().find(|c| c.arc == *ar && c.line == *l && c.horizontal == horizontal).map(|c| c.cid),
333        (ArcCenter(a), ArcCenter(b)) => s.axis_distance_aa_ce_ce.iter().find(|c| c.a == *a && c.b == *b && c.horizontal == horizontal).map(|c| c.cid),
334        (ArcCenter(a), ArcStart(b))  => s.axis_distance_aa_ce_s.iter().find(|c| c.a == *a && c.b == *b && c.horizontal == horizontal).map(|c| c.cid),
335        (ArcCenter(a), ArcEnd(b))    => s.axis_distance_aa_ce_e.iter().find(|c| c.a == *a && c.b == *b && c.horizontal == horizontal).map(|c| c.cid),
336        (ArcStart(a), ArcCenter(b))  => s.axis_distance_aa_s_ce.iter().find(|c| c.a == *a && c.b == *b && c.horizontal == horizontal).map(|c| c.cid),
337        (ArcStart(a), ArcStart(b))   => s.axis_distance_aa_s_s.iter().find(|c| c.a == *a && c.b == *b && c.horizontal == horizontal).map(|c| c.cid),
338        (ArcStart(a), ArcEnd(b))     => s.axis_distance_aa_s_e.iter().find(|c| c.a == *a && c.b == *b && c.horizontal == horizontal).map(|c| c.cid),
339        (ArcEnd(a), ArcCenter(b))    => s.axis_distance_aa_e_ce.iter().find(|c| c.a == *a && c.b == *b && c.horizontal == horizontal).map(|c| c.cid),
340        (ArcEnd(a), ArcStart(b))     => s.axis_distance_aa_e_s.iter().find(|c| c.a == *a && c.b == *b && c.horizontal == horizontal).map(|c| c.cid),
341        (ArcEnd(a), ArcEnd(b))       => s.axis_distance_aa_e_e.iter().find(|c| c.a == *a && c.b == *b && c.horizontal == horizontal).map(|c| c.cid),
342    }
343}
344
345/// Locate the cid of the `PointPointDistance` backing constraint for
346/// a dimension. Each supported endpoint pair maps to one of the
347/// `distance_*` collections (see `arael-sketch/src/actions.rs`
348/// distance dispatch for the canonical routing).
349fn find_distance_cid(
350    s: &Sketch,
351    a: &crate::dimensions::DimensionEndpoint,
352    b: &crate::dimensions::DimensionEndpoint,
353) -> Option<u32> {
354    use crate::dimensions::DimensionEndpoint::*;
355    match (a, b) {
356        (Point(pa), Point(pb)) =>
357            s.distance_pp.iter().find(|c| c.a == *pa && c.b == *pb).map(|c| c.cid),
358        (LineP1(la), LineP1(lb)) => s.distance_ll11.iter().find(|c| c.a == *la && c.b == *lb).map(|c| c.cid),
359        (LineP1(la), LineP2(lb)) => s.distance_ll12.iter().find(|c| c.a == *la && c.b == *lb).map(|c| c.cid),
360        (LineP2(la), LineP1(lb)) => s.distance_ll21.iter().find(|c| c.a == *la && c.b == *lb).map(|c| c.cid),
361        (LineP2(la), LineP2(lb)) => s.distance_ll22.iter().find(|c| c.a == *la && c.b == *lb).map(|c| c.cid),
362        (LineP1(l), Point(p)) | (Point(p), LineP1(l)) =>
363            s.distance_lp1.iter().find(|c| c.line == *l && c.point == *p).map(|c| c.cid),
364        (LineP2(l), Point(p)) | (Point(p), LineP2(l)) =>
365            s.distance_lp2.iter().find(|c| c.line == *l && c.point == *p).map(|c| c.cid),
366        (ArcCenter(ar), Point(p)) | (Point(p), ArcCenter(ar)) =>
367            s.distance_arc_center_p.iter().find(|c| c.arc == *ar && c.point == *p).map(|c| c.cid),
368        (ArcStart(ar), Point(p)) | (Point(p), ArcStart(ar)) =>
369            s.distance_arc_start_p.iter().find(|c| c.arc == *ar && c.point == *p).map(|c| c.cid),
370        (ArcEnd(ar), Point(p)) | (Point(p), ArcEnd(ar)) =>
371            s.distance_arc_end_p.iter().find(|c| c.arc == *ar && c.point == *p).map(|c| c.cid),
372        (ArcCenter(a), ArcCenter(b)) => s.distance_aa_ce_ce.iter().find(|c| c.a == *a && c.b == *b).map(|c| c.cid),
373        (ArcCenter(a), ArcStart(b))  => s.distance_aa_ce_s.iter().find(|c| c.a == *a && c.b == *b).map(|c| c.cid),
374        (ArcCenter(a), ArcEnd(b))    => s.distance_aa_ce_e.iter().find(|c| c.a == *a && c.b == *b).map(|c| c.cid),
375        (ArcStart(a), ArcCenter(b))  => s.distance_aa_s_ce.iter().find(|c| c.a == *a && c.b == *b).map(|c| c.cid),
376        (ArcStart(a), ArcStart(b))   => s.distance_aa_s_s.iter().find(|c| c.a == *a && c.b == *b).map(|c| c.cid),
377        (ArcStart(a), ArcEnd(b))     => s.distance_aa_s_e.iter().find(|c| c.a == *a && c.b == *b).map(|c| c.cid),
378        (ArcEnd(a), ArcCenter(b))    => s.distance_aa_e_ce.iter().find(|c| c.a == *a && c.b == *b).map(|c| c.cid),
379        (ArcEnd(a), ArcStart(b))     => s.distance_aa_e_s.iter().find(|c| c.a == *a && c.b == *b).map(|c| c.cid),
380        (ArcEnd(a), ArcEnd(b))       => s.distance_aa_e_e.iter().find(|c| c.a == *a && c.b == *b).map(|c| c.cid),
381        (LineP1(_), ArcCenter(_)) | (ArcCenter(_), LineP1(_))
382        | (LineP2(_), ArcCenter(_)) | (ArcCenter(_), LineP2(_))
383        | (LineP1(_), ArcStart(_)) | (ArcStart(_), LineP1(_))
384        | (LineP2(_), ArcStart(_)) | (ArcStart(_), LineP2(_))
385        | (LineP1(_), ArcEnd(_)) | (ArcEnd(_), LineP1(_))
386        | (LineP2(_), ArcEnd(_)) | (ArcEnd(_), LineP2(_)) => None,
387    }
388}
389
390/// Locate the cid of the point-to-line distance backing constraint
391/// for a dimension. Covers `distance_pl`, `distance_lp1l`,
392/// `distance_lp2l`, `distance_arc_center_l`, `distance_arc_start_l`,
393/// `distance_arc_end_l`.
394fn find_point_line_cid(
395    s: &Sketch,
396    ep: &crate::dimensions::DimensionEndpoint,
397    line: Ref<Line>,
398) -> Option<u32> {
399    use crate::dimensions::DimensionEndpoint::*;
400    match ep {
401        Point(p) => s.distance_pl.iter().find(|c| c.point == *p && c.line == line).map(|c| c.cid),
402        LineP1(l) => s.distance_lp1l.iter().find(|c| c.a == *l && c.b == line).map(|c| c.cid),
403        LineP2(l) => s.distance_lp2l.iter().find(|c| c.a == *l && c.b == line).map(|c| c.cid),
404        ArcCenter(ar) => s.distance_arc_center_l.iter().find(|c| c.arc == *ar && c.line == line).map(|c| c.cid),
405        ArcStart(ar) => s.distance_arc_start_l.iter().find(|c| c.arc == *ar && c.line == line).map(|c| c.cid),
406        ArcEnd(ar) => s.distance_arc_end_l.iter().find(|c| c.arc == *ar && c.line == line).map(|c| c.cid),
407    }
408}
409
410impl Sketch {
411    /// Create an empty sketch with default solver parameters.
412    pub fn new() -> Self {
413        let drift_sigma = 1000.0_f64;
414        Sketch {
415            points: Arena::new(),
416            lines: Arena::new(),
417            arcs: Arena::new(),
418            drift_isigma: 1.0 / drift_sigma,
419            constraint_isigma: 1000.0, // tight constraints
420            min_length: 0.0001,
421            verbose: false,
422            next_point_id: 0,
423            next_line_id: 0,
424            next_arc_id: 0,
425            coincident_pp: Vec::new(),
426            coincident_lp1: Vec::new(),
427            coincident_lp2: Vec::new(),
428            coincident_ll11: Vec::new(),
429            coincident_ll12: Vec::new(),
430            coincident_ll21: Vec::new(),
431            coincident_ll22: Vec::new(),
432            distance_pp: Vec::new(),
433            hdistance_pp: Vec::new(),
434            vdistance_pp: Vec::new(),
435            point_on_line: Vec::new(),
436            midpoint: Vec::new(),
437            midpoint_lp1: Vec::new(),
438            midpoint_lp2: Vec::new(),
439            midpoint_arc_start: Vec::new(),
440            midpoint_arc_end: Vec::new(),
441            midpoint_arc_point: Vec::new(),
442            midpoint_lp1_arc: Vec::new(),
443            midpoint_lp2_arc: Vec::new(),
444            midpoint_arc_start_arc: Vec::new(),
445            midpoint_arc_end_arc: Vec::new(),
446            point_on_arc: Vec::new(),
447            parallel: Vec::new(),
448            perpendicular: Vec::new(),
449            arc_line_parallel: Vec::new(),
450            arc_arc_parallel: Vec::new(),
451            collinear: Vec::new(),
452            equal_length: Vec::new(),
453            angle: Vec::new(),
454            tangent_la: Vec::new(),
455            concentric: Vec::new(),
456            equal_radius: Vec::new(),
457            tangent_aa: Vec::new(),
458            symmetry_ll: Vec::new(),
459            symmetry_pp: Vec::new(),
460            symmetry_aa: Vec::new(),
461            distance_pl: Vec::new(),
462            distance_lp1l: Vec::new(),
463            distance_lp2l: Vec::new(),
464            distance_arc_center_l: Vec::new(),
465            distance_arc_start_l: Vec::new(),
466            distance_arc_end_l: Vec::new(),
467            line_p1_on_line: Vec::new(),
468            line_p2_on_line: Vec::new(),
469            coincident_arc_center: Vec::new(),
470            coincident_arc_start: Vec::new(),
471            coincident_arc_end: Vec::new(),
472            coincident_lp1_arc_center: Vec::new(),
473            coincident_lp2_arc_center: Vec::new(),
474            coincident_lp1_arc_start: Vec::new(),
475            coincident_lp2_arc_start: Vec::new(),
476            coincident_lp1_arc_end: Vec::new(),
477            coincident_lp2_arc_end: Vec::new(),
478            coincident_arc_center_start: Vec::new(),
479            coincident_arc_center_end: Vec::new(),
480            coincident_arc_start_center: Vec::new(),
481            coincident_arc_end_center: Vec::new(),
482            coincident_arc_start_start: Vec::new(),
483            coincident_arc_start_end: Vec::new(),
484            coincident_arc_end_start: Vec::new(),
485            coincident_arc_end_end: Vec::new(),
486            line_p1_on_arc: Vec::new(),
487            line_p2_on_arc: Vec::new(),
488            distance_ll11: Vec::new(),
489            distance_ll12: Vec::new(),
490            distance_ll21: Vec::new(),
491            distance_ll22: Vec::new(),
492            distance_lp1: Vec::new(),
493            distance_lp2: Vec::new(),
494            distance_arc_center_p: Vec::new(),
495            distance_arc_start_p: Vec::new(),
496            distance_arc_end_p: Vec::new(),
497            distance_arc_center_l1: Vec::new(),
498            distance_arc_center_l2: Vec::new(),
499            distance_arc_start_l1: Vec::new(),
500            distance_arc_start_l2: Vec::new(),
501            distance_arc_end_l1: Vec::new(),
502            distance_arc_end_l2: Vec::new(),
503            distance_aa_ce_ce: Vec::new(),
504            distance_aa_ce_s: Vec::new(),
505            distance_aa_ce_e: Vec::new(),
506            distance_aa_s_ce: Vec::new(),
507            distance_aa_s_s: Vec::new(),
508            distance_aa_s_e: Vec::new(),
509            distance_aa_e_ce: Vec::new(),
510            distance_aa_e_s: Vec::new(),
511            distance_aa_e_e: Vec::new(),
512            distance_concentric: Vec::new(),
513            axis_distance_ll11: Vec::new(),
514            axis_distance_ll12: Vec::new(),
515            axis_distance_ll21: Vec::new(),
516            axis_distance_ll22: Vec::new(),
517            axis_distance_lp1: Vec::new(),
518            axis_distance_lp2: Vec::new(),
519            axis_distance_arc_center_p: Vec::new(),
520            axis_distance_arc_start_p: Vec::new(),
521            axis_distance_arc_end_p: Vec::new(),
522            axis_distance_arc_center_l1: Vec::new(),
523            axis_distance_arc_center_l2: Vec::new(),
524            axis_distance_arc_start_l1: Vec::new(),
525            axis_distance_arc_start_l2: Vec::new(),
526            axis_distance_arc_end_l1: Vec::new(),
527            axis_distance_arc_end_l2: Vec::new(),
528            axis_distance_aa_ce_ce: Vec::new(),
529            axis_distance_aa_ce_s: Vec::new(),
530            axis_distance_aa_ce_e: Vec::new(),
531            axis_distance_aa_s_ce: Vec::new(),
532            axis_distance_aa_s_s: Vec::new(),
533            axis_distance_aa_s_e: Vec::new(),
534            axis_distance_aa_e_ce: Vec::new(),
535            axis_distance_aa_e_s: Vec::new(),
536            axis_distance_aa_e_e: Vec::new(),
537            dimensions: Vec::new(),
538            next_dimension_id: 0,
539            next_constraint_id: 1,
540            user_params: Vec::new(),
541            expr_constraints: Vec::new(),
542            symbol_bag: None,
543            expr_hb: TripletBlock::new(),
544            cached_dof: None,
545        }
546    }
547
548    /// Walk every Vec-stored constraint in canonical order and assign a
549    /// numeric id (C<nid>) to any constraint whose nid is still the 0
550    /// sentinel. Call at the tail of every mutating action and after
551    /// loading a sketch so freshly-deserialised sketches pick up names.
552    pub fn assign_constraint_names(&mut self) {
553        macro_rules! assign {
554            ($($field:ident),* $(,)?) => {
555                $(
556                    for c in self.$field.iter_mut() {
557                        if c.nid == 0 {
558                            c.nid = self.next_constraint_id;
559                            self.next_constraint_id += 1;
560                        }
561                    }
562                )*
563            };
564        }
565        assign!(
566            coincident_pp,
567            coincident_lp1, coincident_lp2,
568            coincident_ll11, coincident_ll12, coincident_ll21, coincident_ll22,
569            distance_pp, hdistance_pp, vdistance_pp,
570            point_on_line,
571            midpoint, midpoint_lp1, midpoint_lp2,
572            midpoint_arc_start, midpoint_arc_end,
573            midpoint_arc_point,
574            midpoint_lp1_arc, midpoint_lp2_arc,
575            midpoint_arc_start_arc, midpoint_arc_end_arc,
576            point_on_arc,
577            parallel, perpendicular, collinear,
578            equal_length, angle,
579            tangent_la, concentric, equal_radius, tangent_aa,
580            symmetry_ll, symmetry_pp, symmetry_aa,
581            distance_pl, distance_lp1l, distance_lp2l,
582            distance_arc_center_l, distance_arc_start_l, distance_arc_end_l,
583            line_p1_on_line, line_p2_on_line,
584            coincident_arc_center, coincident_arc_start, coincident_arc_end,
585            coincident_lp1_arc_center, coincident_lp2_arc_center,
586            coincident_lp1_arc_start, coincident_lp2_arc_start,
587            coincident_lp1_arc_end, coincident_lp2_arc_end,
588            coincident_arc_center_start, coincident_arc_center_end,
589            coincident_arc_start_center, coincident_arc_end_center,
590            coincident_arc_start_start, coincident_arc_start_end,
591            coincident_arc_end_start, coincident_arc_end_end,
592            line_p1_on_arc, line_p2_on_arc,
593            distance_ll11, distance_ll12, distance_ll21, distance_ll22,
594            distance_lp1, distance_lp2,
595            distance_arc_center_p, distance_arc_start_p, distance_arc_end_p,
596            distance_arc_center_l1, distance_arc_center_l2,
597            distance_arc_start_l1, distance_arc_start_l2,
598            distance_arc_end_l1, distance_arc_end_l2,
599            distance_aa_ce_ce, distance_aa_ce_s, distance_aa_ce_e,
600            distance_aa_s_ce, distance_aa_s_s, distance_aa_s_e,
601            distance_aa_e_ce, distance_aa_e_s, distance_aa_e_e,
602            distance_concentric,
603            axis_distance_ll11, axis_distance_ll12, axis_distance_ll21, axis_distance_ll22,
604            axis_distance_lp1, axis_distance_lp2,
605            axis_distance_arc_center_p, axis_distance_arc_start_p, axis_distance_arc_end_p,
606            axis_distance_arc_center_l1, axis_distance_arc_center_l2,
607            axis_distance_arc_start_l1, axis_distance_arc_start_l2,
608            axis_distance_arc_end_l1, axis_distance_arc_end_l2,
609            axis_distance_aa_ce_ce, axis_distance_aa_ce_s, axis_distance_aa_ce_e,
610            axis_distance_aa_s_ce, axis_distance_aa_s_s, axis_distance_aa_s_e,
611            axis_distance_aa_e_ce, axis_distance_aa_e_s, axis_distance_aa_e_e,
612        );
613    }
614
615    /// Add a free point at the given position.
616    pub fn add_point(&mut self, pos: vect2d) -> Ref<Point> {
617        let name = format!("P{}", self.next_point_id);
618        self.next_point_id += 1;
619        self.points.push(Point {
620            pos: Param::new(pos),
621            constraints: PointConstraints { has_fix_x: false, fix_x: 0.0, has_fix_y: false, fix_y: 0.0 },
622            helper: false, quiet: false, name,
623            drag_pull: 0.0,
624            cid: 0, hb: SelfBlock::new(),
625        })
626    }
627
628    /// Add a fixed (non-optimizable) point at the given position.
629    pub fn add_point_fixed(&mut self, pos: vect2d) -> Ref<Point> {
630        let name = format!("P{}", self.next_point_id);
631        self.next_point_id += 1;
632        self.points.push(Point {
633            pos: Param::fixed(pos),
634            constraints: PointConstraints { has_fix_x: false, fix_x: 0.0, has_fix_y: false, fix_y: 0.0 },
635            helper: false, quiet: false, name,
636            drag_pull: 0.0,
637            cid: 0, hb: SelfBlock::new(),
638        })
639    }
640
641    /// Add a helper point (auto-removed when no constraints reference it).
642    pub fn add_helper_point(&mut self, pos: vect2d) -> Ref<Point> {
643        let name = format!("Pc{}", self.next_point_id);
644        self.next_point_id += 1;
645        self.points.push(Point {
646            pos: Param::new(pos),
647            constraints: PointConstraints { has_fix_x: false, fix_x: 0.0, has_fix_y: false, fix_y: 0.0 },
648            helper: true, quiet: false, name,
649            drag_pull: 0.0,
650            cid: 0, hb: SelfBlock::new(),
651        })
652    }
653
654    /// Add a line with two free endpoints.
655    pub fn add_line(&mut self, p1: vect2d, p2: vect2d) -> Ref<Line> {
656        let name = format!("L{}", self.next_line_id);
657        self.next_line_id += 1;
658        self.lines.push(Line {
659            p1: Param::new(p1),
660            p2: Param::new(p2),
661            constraints: LineConstraints { horizontal: false, vertical: false, has_length: false, length: 0.0, has_angle: false, target_angle: 0.0, h_dir_sign: f64::NAN, v_dir_sign: f64::NAN },
662            style: LineStyle::Solid, construction: false, quiet: false, name,
663            cid: 0, hb: SelfBlock::new(),
664        })
665    }
666
667    /// Add an arc or circle. When `closed` is true, start/end angles are
668    /// fixed (not optimized) since they are meaningless for a full circle.
669    pub fn add_arc(&mut self, center: vect2d, radius: f64, start: f64, end: f64, closed: bool) -> Ref<Arc> {
670        self.add_arc_with_dir(center, radius, start, end, closed, true)
671    }
672
673    /// Add an arc with explicit direction. ccw=true means CCW from start to end.
674    /// For CW arcs, end_angle is adjusted so that end - start < 0.
675    pub fn add_arc_with_dir(&mut self, center: vect2d, radius: f64, start: f64, end: f64, closed: bool, ccw: bool) -> Ref<Arc> {
676        // Ensure end - start has the correct sign for the arc direction
677        let end = if !closed && !ccw && end > start {
678            end - std::f64::consts::TAU
679        } else if !closed && ccw && end < start {
680            end + std::f64::consts::TAU
681        } else {
682            end
683        };
684        let name = format!("A{}", self.next_arc_id);
685        self.next_arc_id += 1;
686        self.arcs.push(Arc {
687            center: Param::new(center),
688            radius: Param::new(radius),
689            radius_b: Param::new(radius),
690            rotation: Param::fixed(0.0),
691            start_angle: if closed { Param::fixed(start) } else { Param::new(start) },
692            end_angle: if closed { Param::fixed(end) } else { Param::new(end) },
693            closed, ccw,
694            is_ellipse: false,
695            style: LineStyle::Solid, construction: false, quiet: false, name,
696            constraints: ArcConstraints {
697                has_target_radius: false, target_radius: 0.0,
698                has_target_radius_b: false, target_radius_b: 0.0,
699                has_target_sweep: false, target_sweep: 0.0, sweep_sign: 1.0, has_target_rotation: false, target_rotation: 0.0,
700            },
701            cid: 0, hb: SelfBlock::new(),
702        })
703    }
704
705    /// Add an ellipse (closed) or elliptic arc. rx = semi-major, ry = semi-minor,
706    /// rot = rotation angle of the ellipse axes.
707    pub fn add_ellipse(&mut self, center: vect2d, rx: f64, ry: f64, rot: f64, closed: bool) -> Ref<Arc> {
708        let name = format!("EA{}", self.next_arc_id);
709        self.next_arc_id += 1;
710        self.arcs.push(Arc {
711            center: Param::new(center),
712            radius: Param::new(rx),
713            radius_b: Param::new(ry),
714            rotation: Param::new(rot),
715            start_angle: if closed { Param::fixed(0.0) } else { Param::new(0.0) },
716            end_angle: if closed { Param::fixed(std::f64::consts::TAU) } else { Param::new(std::f64::consts::TAU) },
717            closed,
718            is_ellipse: true,
719            ccw: true,
720            style: LineStyle::Solid, construction: false, quiet: false, name,
721            constraints: ArcConstraints {
722                has_target_radius: false, target_radius: 0.0,
723                has_target_radius_b: false, target_radius_b: 0.0,
724                has_target_sweep: false, target_sweep: 0.0, sweep_sign: 1.0, has_target_rotation: false, target_rotation: 0.0,
725            },
726            cid: 0, hb: SelfBlock::new(),
727        })
728    }
729
730    /// Add a partial elliptic arc with explicit center parameterization.
731    pub fn add_elliptic_arc(&mut self, center: vect2d, rx: f64, ry: f64,
732        rot: f64, start: f64, end: f64, ccw: bool) -> Ref<Arc>
733    {
734        let end = if !ccw && end > start {
735            end - std::f64::consts::TAU
736        } else if ccw && end < start {
737            end + std::f64::consts::TAU
738        } else {
739            end
740        };
741        let name = format!("EA{}", self.next_arc_id);
742        self.next_arc_id += 1;
743        self.arcs.push(Arc {
744            center: Param::new(center),
745            radius: Param::new(rx),
746            radius_b: Param::new(ry),
747            rotation: Param::new(rot),
748            start_angle: Param::new(start),
749            end_angle: Param::new(end),
750            closed: false,
751            is_ellipse: true,
752            ccw,
753            style: LineStyle::Solid, construction: false, quiet: false, name,
754            constraints: ArcConstraints {
755                has_target_radius: false, target_radius: 0.0,
756                has_target_radius_b: false, target_radius_b: 0.0,
757                has_target_sweep: false, target_sweep: 0.0, sweep_sign: 1.0, has_target_rotation: false, target_rotation: 0.0,
758            },
759            cid: 0, hb: SelfBlock::new(),
760        })
761    }
762
763    /// Remove a point and all constraints referencing it.
764    pub fn delete_point(&mut self, r: Ref<Point>) {
765        self.dimensions.retain(|d| !d.kind.references_point(r));
766        self.points.remove(r);
767        self.coincident_pp.retain(|c| c.a != r && c.b != r);
768        self.coincident_lp1.retain(|c| c.point != r);
769        self.coincident_lp2.retain(|c| c.point != r);
770        self.distance_pp.retain(|c| c.a != r && c.b != r);
771        self.hdistance_pp.retain(|c| c.a != r && c.b != r);
772        self.vdistance_pp.retain(|c| c.a != r && c.b != r);
773        self.point_on_line.retain(|c| c.point != r);
774        self.midpoint.retain(|c| c.point != r);
775        self.midpoint_arc_point.retain(|c| c.point != r);
776        self.point_on_arc.retain(|c| c.point != r);
777        self.distance_pl.retain(|c| c.point != r);
778        self.coincident_arc_center.retain(|c| c.point != r);
779        self.coincident_arc_start.retain(|c| c.point != r);
780        self.coincident_arc_end.retain(|c| c.point != r);
781        self.distance_lp1.retain(|c| c.point != r);
782        self.distance_lp2.retain(|c| c.point != r);
783        self.distance_arc_center_p.retain(|c| c.point != r);
784        self.distance_arc_start_p.retain(|c| c.point != r);
785        self.distance_arc_end_p.retain(|c| c.point != r);
786        self.symmetry_pp.retain(|c| c.a != r && c.c != r);
787        self.axis_distance_lp1.retain(|c| c.point != r);
788        self.axis_distance_lp2.retain(|c| c.point != r);
789        self.axis_distance_arc_center_p.retain(|c| c.point != r);
790        self.axis_distance_arc_start_p.retain(|c| c.point != r);
791        self.axis_distance_arc_end_p.retain(|c| c.point != r);
792    }
793
794    /// Remove a line and all constraints referencing it.
795    pub fn delete_line(&mut self, r: Ref<Line>) {
796        self.dimensions.retain(|d| !d.kind.references_line(r));
797        self.lines.remove(r);
798        self.coincident_lp1.retain(|c| c.line != r);
799        self.coincident_lp2.retain(|c| c.line != r);
800        self.coincident_ll11.retain(|c| c.a != r && c.b != r);
801        self.coincident_ll12.retain(|c| c.a != r && c.b != r);
802        self.coincident_ll21.retain(|c| c.a != r && c.b != r);
803        self.coincident_ll22.retain(|c| c.a != r && c.b != r);
804        self.point_on_line.retain(|c| c.line != r);
805        self.midpoint.retain(|c| c.line != r);
806        self.midpoint_lp1.retain(|c| c.line != r && c.target != r);
807        self.midpoint_lp2.retain(|c| c.line != r && c.target != r);
808        self.midpoint_arc_start.retain(|c| c.line != r);
809        self.midpoint_arc_end.retain(|c| c.line != r);
810        self.midpoint_lp1_arc.retain(|c| c.line != r);
811        self.midpoint_lp2_arc.retain(|c| c.line != r);
812        self.parallel.retain(|c| c.a != r && c.b != r);
813        self.perpendicular.retain(|c| c.a != r && c.b != r);
814        self.collinear.retain(|c| c.a != r && c.b != r);
815        self.equal_length.retain(|c| c.a != r && c.b != r);
816        self.angle.retain(|c| c.a != r && c.b != r);
817        self.tangent_la.retain(|c| c.line != r);
818        self.symmetry_ll.retain(|c| c.a != r && c.b != r && c.c != r);
819        self.symmetry_pp.retain(|c| c.line != r);
820        self.symmetry_aa.retain(|c| c.line != r);
821        self.distance_pl.retain(|c| c.line != r);
822        self.distance_lp1l.retain(|c| c.a != r && c.b != r);
823        self.distance_lp2l.retain(|c| c.a != r && c.b != r);
824        self.distance_arc_center_l.retain(|c| c.line != r);
825        self.distance_arc_start_l.retain(|c| c.line != r);
826        self.distance_arc_end_l.retain(|c| c.line != r);
827        self.line_p1_on_line.retain(|c| c.a != r && c.b != r);
828        self.line_p2_on_line.retain(|c| c.a != r && c.b != r);
829        self.coincident_lp1_arc_center.retain(|c| c.line != r);
830        self.coincident_lp2_arc_center.retain(|c| c.line != r);
831        self.coincident_lp1_arc_start.retain(|c| c.line != r);
832        self.coincident_lp2_arc_start.retain(|c| c.line != r);
833        self.coincident_lp1_arc_end.retain(|c| c.line != r);
834        self.coincident_lp2_arc_end.retain(|c| c.line != r);
835        self.line_p1_on_arc.retain(|c| c.line != r);
836        self.line_p2_on_arc.retain(|c| c.line != r);
837        self.distance_ll11.retain(|c| c.a != r && c.b != r);
838        self.distance_ll12.retain(|c| c.a != r && c.b != r);
839        self.distance_ll21.retain(|c| c.a != r && c.b != r);
840        self.distance_ll22.retain(|c| c.a != r && c.b != r);
841        self.distance_lp1.retain(|c| c.line != r);
842        self.distance_lp2.retain(|c| c.line != r);
843        self.distance_arc_center_l1.retain(|c| c.line != r);
844        self.distance_arc_center_l2.retain(|c| c.line != r);
845        self.distance_arc_start_l1.retain(|c| c.line != r);
846        self.distance_arc_start_l2.retain(|c| c.line != r);
847        self.distance_arc_end_l1.retain(|c| c.line != r);
848        self.distance_arc_end_l2.retain(|c| c.line != r);
849        self.axis_distance_ll11.retain(|c| c.a != r && c.b != r);
850        self.axis_distance_ll12.retain(|c| c.a != r && c.b != r);
851        self.axis_distance_ll21.retain(|c| c.a != r && c.b != r);
852        self.axis_distance_ll22.retain(|c| c.a != r && c.b != r);
853        self.axis_distance_lp1.retain(|c| c.line != r);
854        self.axis_distance_lp2.retain(|c| c.line != r);
855        self.axis_distance_arc_center_l1.retain(|c| c.line != r);
856        self.axis_distance_arc_center_l2.retain(|c| c.line != r);
857        self.axis_distance_arc_start_l1.retain(|c| c.line != r);
858        self.axis_distance_arc_start_l2.retain(|c| c.line != r);
859        self.axis_distance_arc_end_l1.retain(|c| c.line != r);
860        self.axis_distance_arc_end_l2.retain(|c| c.line != r);
861        self.cleanup_helper_points();
862    }
863
864    /// Remove an arc and all constraints referencing it.
865    pub fn delete_arc(&mut self, r: Ref<Arc>) {
866        self.dimensions.retain(|d| !d.kind.references_arc(r));
867        self.arcs.remove(r);
868        self.point_on_arc.retain(|c| c.arc != r);
869        self.line_p1_on_arc.retain(|c| c.arc != r);
870        self.line_p2_on_arc.retain(|c| c.arc != r);
871        self.tangent_la.retain(|c| c.arc != r);
872        self.concentric.retain(|c| c.a != r && c.b != r);
873        self.equal_radius.retain(|c| c.a != r && c.b != r);
874        self.tangent_aa.retain(|c| c.a != r && c.b != r);
875        self.symmetry_aa.retain(|c| c.a != r && c.c != r);
876        self.midpoint_arc_start.retain(|c| c.arc != r);
877        self.midpoint_arc_end.retain(|c| c.arc != r);
878        self.midpoint_arc_point.retain(|c| c.arc != r);
879        self.midpoint_lp1_arc.retain(|c| c.arc != r);
880        self.midpoint_lp2_arc.retain(|c| c.arc != r);
881        self.midpoint_arc_start_arc.retain(|c| c.a != r && c.b != r);
882        self.midpoint_arc_end_arc.retain(|c| c.a != r && c.b != r);
883        self.coincident_arc_center.retain(|c| c.arc != r);
884        self.coincident_arc_start.retain(|c| c.arc != r);
885        self.coincident_arc_end.retain(|c| c.arc != r);
886        self.coincident_lp1_arc_center.retain(|c| c.arc != r);
887        self.coincident_lp2_arc_center.retain(|c| c.arc != r);
888        self.coincident_lp1_arc_start.retain(|c| c.arc != r);
889        self.coincident_lp2_arc_start.retain(|c| c.arc != r);
890        self.coincident_lp1_arc_end.retain(|c| c.arc != r);
891        self.coincident_lp2_arc_end.retain(|c| c.arc != r);
892        self.coincident_arc_center_start.retain(|c| c.a != r && c.b != r);
893        self.coincident_arc_center_end.retain(|c| c.a != r && c.b != r);
894        self.coincident_arc_start_center.retain(|c| c.a != r && c.b != r);
895        self.coincident_arc_end_center.retain(|c| c.a != r && c.b != r);
896        self.coincident_arc_start_start.retain(|c| c.a != r && c.b != r);
897        self.coincident_arc_start_end.retain(|c| c.a != r && c.b != r);
898        self.coincident_arc_end_start.retain(|c| c.a != r && c.b != r);
899        self.coincident_arc_end_end.retain(|c| c.a != r && c.b != r);
900        self.distance_arc_center_p.retain(|c| c.arc != r);
901        self.distance_arc_start_p.retain(|c| c.arc != r);
902        self.distance_arc_end_p.retain(|c| c.arc != r);
903        self.distance_arc_center_l.retain(|c| c.arc != r);
904        self.distance_arc_start_l.retain(|c| c.arc != r);
905        self.distance_arc_end_l.retain(|c| c.arc != r);
906        self.distance_arc_center_l1.retain(|c| c.arc != r);
907        self.distance_arc_center_l2.retain(|c| c.arc != r);
908        self.distance_arc_start_l1.retain(|c| c.arc != r);
909        self.distance_arc_start_l2.retain(|c| c.arc != r);
910        self.distance_arc_end_l1.retain(|c| c.arc != r);
911        self.distance_arc_end_l2.retain(|c| c.arc != r);
912        self.distance_aa_ce_ce.retain(|c| c.a != r && c.b != r);
913        self.distance_aa_ce_s.retain(|c| c.a != r && c.b != r);
914        self.distance_aa_ce_e.retain(|c| c.a != r && c.b != r);
915        self.distance_aa_s_ce.retain(|c| c.a != r && c.b != r);
916        self.distance_aa_s_s.retain(|c| c.a != r && c.b != r);
917        self.distance_aa_s_e.retain(|c| c.a != r && c.b != r);
918        self.distance_aa_e_ce.retain(|c| c.a != r && c.b != r);
919        self.distance_aa_e_s.retain(|c| c.a != r && c.b != r);
920        self.distance_aa_e_e.retain(|c| c.a != r && c.b != r);
921        self.distance_concentric.retain(|c| c.a != r && c.b != r);
922        self.axis_distance_arc_center_p.retain(|c| c.arc != r);
923        self.axis_distance_arc_start_p.retain(|c| c.arc != r);
924        self.axis_distance_arc_end_p.retain(|c| c.arc != r);
925        self.axis_distance_arc_center_l1.retain(|c| c.arc != r);
926        self.axis_distance_arc_center_l2.retain(|c| c.arc != r);
927        self.axis_distance_arc_start_l1.retain(|c| c.arc != r);
928        self.axis_distance_arc_start_l2.retain(|c| c.arc != r);
929        self.axis_distance_arc_end_l1.retain(|c| c.arc != r);
930        self.axis_distance_arc_end_l2.retain(|c| c.arc != r);
931        self.axis_distance_aa_ce_ce.retain(|c| c.a != r && c.b != r);
932        self.axis_distance_aa_ce_s.retain(|c| c.a != r && c.b != r);
933        self.axis_distance_aa_ce_e.retain(|c| c.a != r && c.b != r);
934        self.axis_distance_aa_s_ce.retain(|c| c.a != r && c.b != r);
935        self.axis_distance_aa_s_s.retain(|c| c.a != r && c.b != r);
936        self.axis_distance_aa_s_e.retain(|c| c.a != r && c.b != r);
937        self.axis_distance_aa_e_ce.retain(|c| c.a != r && c.b != r);
938        self.axis_distance_aa_e_s.retain(|c| c.a != r && c.b != r);
939        self.axis_distance_aa_e_e.retain(|c| c.a != r && c.b != r);
940        self.cleanup_helper_points();
941    }
942
943    /// Remove helper points that are no longer needed.
944    /// A helper is removed if it lost its bridge constraint (semantic origin
945    /// gone) or has no purpose constraint. Cascades until stable.
946    pub fn cleanup_helper_points(&mut self) {
947        loop {
948            // Find which helpers have a bridge (know what they represent)
949            let mut has_bridge: std::collections::HashSet<u32> = std::collections::HashSet::new();
950            for c in &self.coincident_lp1 { if self.points.get(c.point).is_some_and(|p| p.helper) { has_bridge.insert(c.point.index()); } }
951            for c in &self.coincident_lp2 { if self.points.get(c.point).is_some_and(|p| p.helper) { has_bridge.insert(c.point.index()); } }
952            for c in &self.coincident_arc_center { if self.points.get(c.point).is_some_and(|p| p.helper) { has_bridge.insert(c.point.index()); } }
953            for c in &self.coincident_arc_start { if self.points.get(c.point).is_some_and(|p| p.helper) { has_bridge.insert(c.point.index()); } }
954            for c in &self.coincident_arc_end { if self.points.get(c.point).is_some_and(|p| p.helper) { has_bridge.insert(c.point.index()); } }
955            for c in &self.coincident_pp {
956                if self.points.get(c.a).is_some_and(|p| p.helper) { has_bridge.insert(c.a.index()); }
957                if self.points.get(c.b).is_some_and(|p| p.helper) { has_bridge.insert(c.b.index()); }
958            }
959
960            // Find which helpers have a purpose constraint
961            let mut has_purpose: std::collections::HashSet<u32> = std::collections::HashSet::new();
962            let mut add_pt = |r: Ref<Point>| { has_purpose.insert(r.index()); };
963            for c in &self.distance_pp { add_pt(c.a); add_pt(c.b); }
964            for c in &self.hdistance_pp { add_pt(c.a); add_pt(c.b); }
965            for c in &self.vdistance_pp { add_pt(c.a); add_pt(c.b); }
966            for c in &self.point_on_line { add_pt(c.point); }
967            for c in &self.midpoint { add_pt(c.point); }
968            for c in &self.point_on_arc { add_pt(c.point); }
969            for c in &self.distance_pl { add_pt(c.point); }
970            for c in &self.distance_lp1 { add_pt(c.point); }
971            for c in &self.distance_lp2 { add_pt(c.point); }
972            for c in &self.symmetry_pp { add_pt(c.a); add_pt(c.c); }
973
974            // Remove helpers that lost their bridge OR have no purpose
975            let to_remove: std::vec::Vec<Ref<Point>> = self.points.refs()
976                .filter(|r| self.points[*r].helper
977                    && (!has_bridge.contains(&r.index()) || !has_purpose.contains(&r.index())))
978                .collect();
979            if to_remove.is_empty() { break; }
980
981            for r in &to_remove {
982                self.coincident_pp.retain(|c| c.a != *r && c.b != *r);
983                self.coincident_lp1.retain(|c| c.point != *r);
984                self.coincident_lp2.retain(|c| c.point != *r);
985                self.coincident_arc_center.retain(|c| c.point != *r);
986                self.coincident_arc_start.retain(|c| c.point != *r);
987                self.coincident_arc_end.retain(|c| c.point != *r);
988                self.symmetry_pp.retain(|c| c.a != *r && c.c != *r);
989                self.distance_pp.retain(|c| c.a != *r && c.b != *r);
990                self.hdistance_pp.retain(|c| c.a != *r && c.b != *r);
991                self.vdistance_pp.retain(|c| c.a != *r && c.b != *r);
992                self.point_on_line.retain(|c| c.point != *r);
993                self.midpoint.retain(|c| c.point != *r);
994                self.point_on_arc.retain(|c| c.point != *r);
995                self.distance_pl.retain(|c| c.point != *r);
996                self.distance_lp1.retain(|c| c.point != *r);
997                self.distance_lp2.retain(|c| c.point != *r);
998            }
999            for r in to_remove { self.points.remove(r); }
1000        }
1001    }
1002
1003    /// Remove duplicate constraints from all collections. Prints a warning if any are found.
1004    pub fn dedup_constraints(&mut self) {
1005        let mut total_removed = 0usize;
1006        macro_rules! dedup_ab {
1007            ($coll:expr, $name:expr, $container_a:expr, $container_b:expr) => {
1008                let mut seen = std::collections::HashSet::new();
1009                for c in $coll.iter() {
1010                    if !seen.insert((c.a.index(), c.b.index())) {
1011                        let na = $container_a.get(c.a).map(|e| e.name.as_str()).unwrap_or("?");
1012                        let nb = $container_b.get(c.b).map(|e| e.name.as_str()).unwrap_or("?");
1013                        eprintln!("BUG: duplicate {} constraint: a={}, b={}", $name, na, nb);
1014                        eprintln!("{}", std::backtrace::Backtrace::force_capture());
1015                        total_removed += 1;
1016                    }
1017                }
1018                seen.clear();
1019                $coll.retain(|c| seen.insert((c.a.index(), c.b.index())));
1020            };
1021        }
1022        macro_rules! dedup_lp {
1023            ($coll:expr, $name:expr) => {
1024                let mut seen = std::collections::HashSet::new();
1025                for c in $coll.iter() {
1026                    if !seen.insert((c.line.index(), c.point.index())) {
1027                        let nl = self.lines.get(c.line).map(|e| e.name.as_str()).unwrap_or("?");
1028                        let np = self.points.get(c.point).map(|e| e.name.as_str()).unwrap_or("?");
1029                        eprintln!("BUG: duplicate {} constraint: line={}, point={}", $name, nl, np);
1030                        eprintln!("{}", std::backtrace::Backtrace::force_capture());
1031                        total_removed += 1;
1032                    }
1033                }
1034                seen.clear();
1035                $coll.retain(|c| seen.insert((c.line.index(), c.point.index())));
1036            };
1037        }
1038        macro_rules! dedup_la {
1039            ($coll:expr, $name:expr) => {
1040                let mut seen = std::collections::HashSet::new();
1041                for c in $coll.iter() {
1042                    if !seen.insert((c.line.index(), c.arc.index())) {
1043                        let nl = self.lines.get(c.line).map(|e| e.name.as_str()).unwrap_or("?");
1044                        let na = self.arcs.get(c.arc).map(|e| e.name.as_str()).unwrap_or("?");
1045                        eprintln!("BUG: duplicate {} constraint: line={}, arc={}", $name, nl, na);
1046                        eprintln!("{}", std::backtrace::Backtrace::force_capture());
1047                        total_removed += 1;
1048                    }
1049                }
1050                seen.clear();
1051                $coll.retain(|c| seen.insert((c.line.index(), c.arc.index())));
1052            };
1053        }
1054        macro_rules! dedup_pa {
1055            ($coll:expr, $name:expr) => {
1056                let mut seen = std::collections::HashSet::new();
1057                for c in $coll.iter() {
1058                    if !seen.insert((c.point.index(), c.arc.index())) {
1059                        let np = self.points.get(c.point).map(|e| e.name.as_str()).unwrap_or("?");
1060                        let na = self.arcs.get(c.arc).map(|e| e.name.as_str()).unwrap_or("?");
1061                        eprintln!("BUG: duplicate {} constraint: point={}, arc={}", $name, np, na);
1062                        eprintln!("{}", std::backtrace::Backtrace::force_capture());
1063                        total_removed += 1;
1064                    }
1065                }
1066                seen.clear();
1067                $coll.retain(|c| seen.insert((c.point.index(), c.arc.index())));
1068            };
1069        }
1070        macro_rules! dedup_pl {
1071            ($coll:expr, $name:expr) => {
1072                let mut seen = std::collections::HashSet::new();
1073                for c in $coll.iter() {
1074                    if !seen.insert((c.point.index(), c.line.index())) {
1075                        let np = self.points.get(c.point).map(|e| e.name.as_str()).unwrap_or("?");
1076                        let nl = self.lines.get(c.line).map(|e| e.name.as_str()).unwrap_or("?");
1077                        eprintln!("BUG: duplicate {} constraint: point={}, line={}", $name, np, nl);
1078                        eprintln!("{}", std::backtrace::Backtrace::force_capture());
1079                        total_removed += 1;
1080                    }
1081                }
1082                seen.clear();
1083                $coll.retain(|c| seen.insert((c.point.index(), c.line.index())));
1084            };
1085        }
1086        dedup_ab!(self.coincident_pp, "coincident_pp", self.points, self.points);
1087        // PP is symmetric: (a,b) == (b,a)
1088        {
1089            let before = self.coincident_pp.len();
1090            let mut seen = std::collections::HashSet::new();
1091            self.coincident_pp.retain(|c| {
1092                let (a, b) = (c.a.index().min(c.b.index()), c.a.index().max(c.b.index()));
1093                seen.insert((a, b))
1094            });
1095            let removed = before - self.coincident_pp.len();
1096            if removed > 0 { eprintln!("BUG: removed {} cross-duplicate coincident_pp constraints", removed); total_removed += removed; }
1097        }
1098        dedup_lp!(self.coincident_lp1, "coincident_lp1");
1099        dedup_lp!(self.coincident_lp2, "coincident_lp2");
1100        dedup_ab!(self.coincident_ll11, "coincident_ll11", self.lines, self.lines);
1101        dedup_ab!(self.coincident_ll12, "coincident_ll12", self.lines, self.lines);
1102        dedup_ab!(self.coincident_ll21, "coincident_ll21", self.lines, self.lines);
1103        dedup_ab!(self.coincident_ll22, "coincident_ll22", self.lines, self.lines);
1104        // Cross-Vec dedup for LL: ll11(a,b)==ll11(b,a), ll22(a,b)==ll22(b,a), ll12(a,b)==ll21(b,a)
1105        {
1106            let mut seen = std::collections::HashSet::new();
1107            // Normalize: represent each endpoint pair as (min_id, max_id) where id encodes line+endpoint
1108            let ep_id = |line: u32, is_p2: bool| -> u64 { (line as u64) << 1 | (is_p2 as u64) };
1109            let mut add = |line_a: u32, p2_a: bool, line_b: u32, p2_b: bool| -> bool {
1110                let a = ep_id(line_a, p2_a);
1111                let b = ep_id(line_b, p2_b);
1112                let key = (a.min(b), a.max(b));
1113                seen.insert(key)
1114            };
1115            let before = self.coincident_ll11.len() + self.coincident_ll12.len()
1116                + self.coincident_ll21.len() + self.coincident_ll22.len();
1117            self.coincident_ll11.retain(|c| add(c.a.index(), false, c.b.index(), false));
1118            self.coincident_ll12.retain(|c| add(c.a.index(), false, c.b.index(), true));
1119            self.coincident_ll21.retain(|c| add(c.a.index(), true, c.b.index(), false));
1120            self.coincident_ll22.retain(|c| add(c.a.index(), true, c.b.index(), true));
1121            let after = self.coincident_ll11.len() + self.coincident_ll12.len()
1122                + self.coincident_ll21.len() + self.coincident_ll22.len();
1123            let removed = before - after;
1124            if removed > 0 { eprintln!("BUG: removed {} cross-duplicate LL coincident constraints", removed); total_removed += removed; }
1125        }
1126        dedup_ab!(self.distance_pp, "distance_pp", self.points, self.points);
1127        dedup_ab!(self.hdistance_pp, "hdistance_pp", self.points, self.points);
1128        dedup_ab!(self.vdistance_pp, "vdistance_pp", self.points, self.points);
1129        dedup_pl!(self.point_on_line, "point_on_line");
1130        dedup_pl!(self.midpoint, "midpoint");
1131        {
1132            let mut seen = std::collections::HashSet::new();
1133            self.midpoint_lp1.retain(|c| seen.insert((c.line.index(), c.target.index())));
1134        }
1135        {
1136            let mut seen = std::collections::HashSet::new();
1137            self.midpoint_lp2.retain(|c| seen.insert((c.line.index(), c.target.index())));
1138        }
1139        {
1140            let mut seen = std::collections::HashSet::new();
1141            self.midpoint_arc_start.retain(|c| seen.insert((c.arc.index(), c.line.index())));
1142        }
1143        {
1144            let mut seen = std::collections::HashSet::new();
1145            self.midpoint_arc_end.retain(|c| seen.insert((c.arc.index(), c.line.index())));
1146        }
1147        {
1148            let mut seen = std::collections::HashSet::new();
1149            self.midpoint_arc_point.retain(|c| seen.insert((c.point.index(), c.arc.index())));
1150        }
1151        {
1152            let mut seen = std::collections::HashSet::new();
1153            self.midpoint_lp1_arc.retain(|c| seen.insert((c.line.index(), c.arc.index())));
1154        }
1155        {
1156            let mut seen = std::collections::HashSet::new();
1157            self.midpoint_lp2_arc.retain(|c| seen.insert((c.line.index(), c.arc.index())));
1158        }
1159        {
1160            let mut seen = std::collections::HashSet::new();
1161            self.midpoint_arc_start_arc.retain(|c| seen.insert((c.a.index(), c.b.index())));
1162        }
1163        {
1164            let mut seen = std::collections::HashSet::new();
1165            self.midpoint_arc_end_arc.retain(|c| seen.insert((c.a.index(), c.b.index())));
1166        }
1167        dedup_pa!(self.point_on_arc, "point_on_arc");
1168        dedup_ab!(self.parallel, "parallel", self.lines, self.lines);
1169        dedup_la!(self.arc_line_parallel, "arc_line_parallel");
1170        dedup_ab!(self.arc_arc_parallel, "arc_arc_parallel", self.arcs, self.arcs);
1171        dedup_ab!(self.perpendicular, "perpendicular", self.lines, self.lines);
1172        dedup_ab!(self.collinear, "collinear", self.lines, self.lines);
1173        {
1174            let mut seen = std::collections::HashSet::new();
1175            self.symmetry_ll.retain(|c| seen.insert((c.a.index(), c.b.index(), c.c.index())));
1176        }
1177        dedup_ab!(self.equal_length, "equal_length", self.lines, self.lines);
1178        dedup_ab!(self.angle, "angle", self.lines, self.lines);
1179        dedup_la!(self.tangent_la, "tangent_la");
1180        dedup_la!(self.line_p1_on_arc, "line_p1_on_arc");
1181        dedup_la!(self.line_p2_on_arc, "line_p2_on_arc");
1182        dedup_ab!(self.concentric, "concentric", self.arcs, self.arcs);
1183        dedup_ab!(self.equal_radius, "equal_radius", self.arcs, self.arcs);
1184        dedup_ab!(self.tangent_aa, "tangent_aa", self.arcs, self.arcs);
1185        dedup_pl!(self.distance_pl, "distance_pl");
1186        dedup_ab!(self.line_p1_on_line, "line_p1_on_line", self.lines, self.lines);
1187        dedup_ab!(self.line_p2_on_line, "line_p2_on_line", self.lines, self.lines);
1188        dedup_pa!(self.coincident_arc_center, "coincident_arc_center");
1189        dedup_pa!(self.coincident_arc_start, "coincident_arc_start");
1190        dedup_pa!(self.coincident_arc_end, "coincident_arc_end");
1191        dedup_la!(self.coincident_lp1_arc_center, "coincident_lp1_arc_center");
1192        dedup_la!(self.coincident_lp2_arc_center, "coincident_lp2_arc_center");
1193        dedup_la!(self.coincident_lp1_arc_start, "coincident_lp1_arc_start");
1194        dedup_la!(self.coincident_lp2_arc_start, "coincident_lp2_arc_start");
1195        dedup_la!(self.coincident_lp1_arc_end, "coincident_lp1_arc_end");
1196        dedup_la!(self.coincident_lp2_arc_end, "coincident_lp2_arc_end");
1197        dedup_ab!(self.coincident_arc_center_start, "coincident_arc_center_start", self.arcs, self.arcs);
1198        dedup_ab!(self.coincident_arc_center_end, "coincident_arc_center_end", self.arcs, self.arcs);
1199        dedup_ab!(self.coincident_arc_start_center, "coincident_arc_start_center", self.arcs, self.arcs);
1200        dedup_ab!(self.coincident_arc_end_center, "coincident_arc_end_center", self.arcs, self.arcs);
1201        dedup_ab!(self.coincident_arc_start_start, "coincident_arc_start_start", self.arcs, self.arcs);
1202        dedup_ab!(self.coincident_arc_start_end, "coincident_arc_start_end", self.arcs, self.arcs);
1203        dedup_ab!(self.coincident_arc_end_start, "coincident_arc_end_start", self.arcs, self.arcs);
1204        dedup_ab!(self.coincident_arc_end_end, "coincident_arc_end_end", self.arcs, self.arcs);
1205        let _ = total_removed;
1206    }
1207
1208    /// Recompute tangent_la sign fields from current geometry.
1209    /// Needed after loading old saves that default sign to 1.0.
1210    pub fn fixup_tangent_signs(&mut self) {
1211        for t in &mut self.tangent_la {
1212            let l = &self.lines[t.line];
1213            let a = &self.arcs[t.arc];
1214            let dx = l.p2.value.x - l.p1.value.x;
1215            let dy = l.p2.value.y - l.p1.value.y;
1216            let len = (dx * dx + dy * dy).sqrt();
1217            if len < 1e-12 { continue; }
1218            let dist = ((a.center.value.x - l.p1.value.x) * dy
1219                      - (a.center.value.y - l.p1.value.y) * dx) / len;
1220            t.sign = if dist >= 0.0 { 1.0 } else { -1.0 };
1221        }
1222    }
1223
1224    /// Merge duplicate helper points at the same position and consolidate
1225    /// helper-point-bridged constraints into direct constraints.
1226    pub fn consolidate_helper_constraints(&mut self) {
1227        // Phase 1: Merge duplicate helper points at the same position.
1228        // If two helper points are at the same position, rewrite all constraints
1229        // referencing the second to reference the first, then remove the second.
1230        let helper_refs: std::vec::Vec<Ref<Point>> = self.points.refs()
1231            .filter(|r| self.points[*r].helper)
1232            .collect();
1233        let mut merged = std::collections::HashMap::<u32, Ref<Point>>::new(); // old -> canonical
1234        for i in 0..helper_refs.len() {
1235            let ri = helper_refs[i];
1236            if merged.contains_key(&ri.index()) { continue; }
1237            let pi = self.points[ri].pos.value;
1238            for j in (i+1)..helper_refs.len() {
1239                let rj = helper_refs[j];
1240                if merged.contains_key(&rj.index()) { continue; }
1241                let pj = self.points[rj].pos.value;
1242                if (pi.x - pj.x).abs() < 1e-9 && (pi.y - pj.y).abs() < 1e-9 {
1243                    merged.insert(rj.index(), ri);
1244                    eprintln!("INFO: merging duplicate helper point {} into {}", rj.index(), ri.index());
1245                }
1246            }
1247        }
1248        if !merged.is_empty() {
1249            // Rewrite all point refs in constraints
1250            let remap = |r: &mut Ref<Point>| {
1251                if let Some(canonical) = merged.get(&r.index()) { *r = *canonical; }
1252            };
1253            for c in &mut self.coincident_pp { remap(&mut c.a); remap(&mut c.b); }
1254            for c in &mut self.coincident_lp1 { remap(&mut c.point); }
1255            for c in &mut self.coincident_lp2 { remap(&mut c.point); }
1256            for c in &mut self.distance_pp { remap(&mut c.a); remap(&mut c.b); }
1257            for c in &mut self.hdistance_pp { remap(&mut c.a); remap(&mut c.b); }
1258            for c in &mut self.vdistance_pp { remap(&mut c.a); remap(&mut c.b); }
1259            for c in &mut self.point_on_line { remap(&mut c.point); }
1260            for c in &mut self.midpoint { remap(&mut c.point); }
1261            for c in &mut self.point_on_arc { remap(&mut c.point); }
1262            for c in &mut self.distance_pl { remap(&mut c.point); }
1263            for c in &mut self.coincident_arc_center { remap(&mut c.point); }
1264            for c in &mut self.coincident_arc_start { remap(&mut c.point); }
1265            for c in &mut self.coincident_arc_end { remap(&mut c.point); }
1266            for c in &mut self.distance_arc_center_p { remap(&mut c.point); }
1267            for c in &mut self.distance_arc_start_p { remap(&mut c.point); }
1268            for c in &mut self.distance_arc_end_p { remap(&mut c.point); }
1269            for c in &mut self.axis_distance_lp1 { remap(&mut c.point); }
1270            for c in &mut self.axis_distance_lp2 { remap(&mut c.point); }
1271            for c in &mut self.axis_distance_arc_center_p { remap(&mut c.point); }
1272            for c in &mut self.axis_distance_arc_start_p { remap(&mut c.point); }
1273            for c in &mut self.axis_distance_arc_end_p { remap(&mut c.point); }
1274            for c in &mut self.symmetry_pp { remap(&mut c.a); remap(&mut c.c); }
1275            for c in &mut self.distance_lp1 { remap(&mut c.point); }
1276            for c in &mut self.distance_lp2 { remap(&mut c.point); }
1277            // Remove merged points
1278            for old in merged.keys() { self.points.remove(Ref::new(*old)); }
1279            // Dedup again after remapping
1280            self.dedup_constraints();
1281        }
1282
1283        // Phase 2: Replace helper-point bridges with direct constraints.
1284        // Collect ALL matching arcs (not just the first) since point merging
1285        // in Phase 1 can create multiple arc constraints per helper point.
1286        let helper_refs: std::vec::Vec<Ref<Point>> = self.points.refs()
1287            .filter(|r| self.points[*r].helper)
1288            .collect();
1289        for hr in &helper_refs {
1290            let hr = *hr;
1291            let lp1s: std::vec::Vec<Ref<Line>> = self.coincident_lp1.iter().filter(|c| c.point == hr).map(|c| c.line).collect();
1292            let lp2s: std::vec::Vec<Ref<Line>> = self.coincident_lp2.iter().filter(|c| c.point == hr).map(|c| c.line).collect();
1293            let acs: std::vec::Vec<Ref<Arc>> = self.coincident_arc_center.iter().filter(|c| c.point == hr).map(|c| c.arc).collect();
1294            let a_starts: std::vec::Vec<Ref<Arc>> = self.coincident_arc_start.iter().filter(|c| c.point == hr).map(|c| c.arc).collect();
1295            let a_ends: std::vec::Vec<Ref<Arc>> = self.coincident_arc_end.iter().filter(|c| c.point == hr).map(|c| c.arc).collect();
1296
1297            macro_rules! consolidate {
1298                ($lines:expr, $arcs:expr, $lp_coll:ident, $arc_coll:ident, $direct_coll:ident, $DirectType:ident, $label:expr) => {
1299                    for &line in &$lines {
1300                        for &arc in &$arcs {
1301                            self.$direct_coll.push($DirectType { line, arc, nid: 0, cid: 0, hb: CrossBlock::new() });
1302                            self.$lp_coll.retain(|c| !(c.line == line && c.point == hr));
1303                            self.$arc_coll.retain(|c| !(c.point == hr && c.arc == arc));
1304                            eprintln!("INFO: consolidated helper {} -> {}", hr.index(), $label);
1305                        }
1306                    }
1307                };
1308            }
1309            consolidate!(lp1s, acs, coincident_lp1, coincident_arc_center, coincident_lp1_arc_center, CoincidentLP1ArcCenter, "LP1ArcCenter");
1310            consolidate!(lp2s, acs, coincident_lp2, coincident_arc_center, coincident_lp2_arc_center, CoincidentLP2ArcCenter, "LP2ArcCenter");
1311            consolidate!(lp1s, a_starts, coincident_lp1, coincident_arc_start, coincident_lp1_arc_start, CoincidentLP1ArcStart, "LP1ArcStart");
1312            consolidate!(lp2s, a_starts, coincident_lp2, coincident_arc_start, coincident_lp2_arc_start, CoincidentLP2ArcStart, "LP2ArcStart");
1313            consolidate!(lp1s, a_ends, coincident_lp1, coincident_arc_end, coincident_lp1_arc_end, CoincidentLP1ArcEnd, "LP1ArcEnd");
1314            consolidate!(lp2s, a_ends, coincident_lp2, coincident_arc_end, coincident_lp2_arc_end, CoincidentLP2ArcEnd, "LP2ArcEnd");
1315        }
1316        self.cleanup_helper_points();
1317        self.dedup_constraints();
1318    }
1319
1320}
1321
1322impl arael::model::ExtendedModel for Sketch {
1323    fn extended_deserialize64(&mut self) {
1324        // After solver writes back optimized values, sync radius_b.value for non-ellipses
1325        let refs: Vec<_> = self.arcs.refs().collect();
1326        for r in refs {
1327            if !self.arcs[r].is_ellipse {
1328                self.arcs[r].radius_b.value = self.arcs[r].radius.value;
1329            }
1330        }
1331    }
1332
1333    fn extended_update64(&mut self, _params: &[f64]) {
1334        // For non-ellipse arcs, keep radius_b work value in sync with radius
1335        let refs: Vec<_> = self.arcs.refs().collect();
1336        for r in refs {
1337            if !self.arcs[r].is_ellipse {
1338                let rv = self.arcs[r].radius.work();
1339                *self.arcs[r].radius_b.work_mut() = rv;
1340            }
1341        }
1342    }
1343
1344    fn extended_cost64(&self, params: &[f64]) -> f64 {
1345        if self.expr_constraints.is_empty() { return 0.0; }
1346        let bag = self.symbol_bag.as_ref().expect("symbol_bag not built");
1347        let vars = bag.eval_vars(params);
1348        let isigma = self.constraint_isigma;
1349        let mut total = 0.0;
1350        for ec in &self.expr_constraints {
1351            match ec.cost(&vars, isigma) {
1352                Ok(c) => total += c,
1353                Err(e) => eprintln!("expr constraint eval error: {}: {}", ec.description, e),
1354            }
1355        }
1356        total
1357    }
1358
1359    fn extended_compute64(&mut self, params: &[f64], grad: &mut [f64]) {
1360        if self.expr_constraints.is_empty() { return; }
1361        let bag = self.symbol_bag.as_ref().expect("symbol_bag not built");
1362        let vars = bag.eval_vars(params);
1363        let isigma = self.constraint_isigma;
1364        let hb = &mut self.expr_hb as *mut TripletBlock<f64>;
1365        for ec in &self.expr_constraints {
1366            if let Err(e) = ec.compute(&vars, isigma, unsafe { &mut *hb }, grad) {
1367                eprintln!("expr constraint eval error: {}: {}", ec.description, e);
1368            }
1369        }
1370    }
1371
1372    fn extended_jacobian64(&mut self, params: &[f64], rows: &mut std::vec::Vec<arael::model::JacobianRow<f64>>, cid: &mut u32) {
1373        if self.expr_constraints.is_empty() { return; }
1374        let bag = self.symbol_bag.as_ref().expect("symbol_bag not built");
1375        let vars = bag.eval_vars(params);
1376        let isigma = self.constraint_isigma;
1377        for ec in &mut self.expr_constraints {
1378            ec.cid = *cid;
1379            match ec.jacobian_row(&vars, isigma) {
1380                Ok((residual, entries)) => {
1381                    rows.push(arael::model::JacobianRow { constraint: *cid, label: ec.label, residual, entries });
1382                }
1383                Err(e) => eprintln!("expr constraint eval error: {}: {}", ec.description, e),
1384            }
1385            *cid += 1;
1386        }
1387    }
1388}
1389
1390impl Sketch {
1391    /// Return `(nid, cid)` pairs for every constraint in the sketch
1392    /// (across all constraint collections, entity constraints excluded).
1393    ///
1394    /// `nid` is the user-visible stable numeric constraint id (C1, C2,
1395    /// ...) and is serialised, so pre/post action comparisons are
1396    /// robust against the macro's per-pass cid renumbering. `cid` is
1397    /// the post-calc_jacobian identifier used by [`Jacobian::rows`].
1398    ///
1399    /// Must be called AFTER `calc_jacobian`/`solve`/`compute_dof` so
1400    /// the `cid` field is populated on every constraint instance.
1401    /// Constraints whose `nid` is still 0 (never named) are omitted.
1402    pub fn constraint_nid_cid_pairs(&self) -> std::vec::Vec<(u32, u32)> {
1403        let mut out = std::vec::Vec::new();
1404        macro_rules! collect {
1405            ($($field:ident),* $(,)?) => {
1406                $(
1407                    for c in &self.$field {
1408                        if c.nid != 0 {
1409                            out.push((c.nid, c.cid));
1410                        }
1411                    }
1412                )*
1413            };
1414        }
1415        collect!(
1416            coincident_pp,
1417            coincident_lp1, coincident_lp2,
1418            coincident_ll11, coincident_ll12, coincident_ll21, coincident_ll22,
1419            distance_pp, hdistance_pp, vdistance_pp,
1420            point_on_line,
1421            midpoint, midpoint_lp1, midpoint_lp2,
1422            midpoint_arc_start, midpoint_arc_end,
1423            midpoint_arc_point,
1424            midpoint_lp1_arc, midpoint_lp2_arc,
1425            midpoint_arc_start_arc, midpoint_arc_end_arc,
1426            point_on_arc,
1427            parallel, perpendicular, collinear,
1428            equal_length, angle,
1429            tangent_la, concentric, equal_radius, tangent_aa,
1430            symmetry_ll, symmetry_pp, symmetry_aa,
1431            distance_pl, distance_lp1l, distance_lp2l,
1432            distance_arc_center_l, distance_arc_start_l, distance_arc_end_l,
1433            line_p1_on_line, line_p2_on_line,
1434            coincident_arc_center, coincident_arc_start, coincident_arc_end,
1435            coincident_lp1_arc_center, coincident_lp2_arc_center,
1436            coincident_lp1_arc_start, coincident_lp2_arc_start,
1437            coincident_lp1_arc_end, coincident_lp2_arc_end,
1438            coincident_arc_center_start, coincident_arc_center_end,
1439            coincident_arc_start_center, coincident_arc_end_center,
1440            coincident_arc_start_start, coincident_arc_start_end,
1441            coincident_arc_end_start, coincident_arc_end_end,
1442            line_p1_on_arc, line_p2_on_arc,
1443            distance_ll11, distance_ll12, distance_ll21, distance_ll22,
1444            distance_lp1, distance_lp2,
1445            distance_arc_center_p, distance_arc_start_p, distance_arc_end_p,
1446            distance_arc_center_l1, distance_arc_center_l2,
1447            distance_arc_start_l1, distance_arc_start_l2,
1448            distance_arc_end_l1, distance_arc_end_l2,
1449            distance_aa_ce_ce, distance_aa_ce_s, distance_aa_ce_e,
1450            distance_aa_s_ce, distance_aa_s_s, distance_aa_s_e,
1451            distance_aa_e_ce, distance_aa_e_s, distance_aa_e_e,
1452            distance_concentric,
1453            axis_distance_ll11, axis_distance_ll12, axis_distance_ll21, axis_distance_ll22,
1454            axis_distance_lp1, axis_distance_lp2,
1455            axis_distance_arc_center_p, axis_distance_arc_start_p, axis_distance_arc_end_p,
1456            axis_distance_arc_center_l1, axis_distance_arc_center_l2,
1457            axis_distance_arc_start_l1, axis_distance_arc_start_l2,
1458            axis_distance_arc_end_l1, axis_distance_arc_end_l2,
1459            axis_distance_aa_ce_ce, axis_distance_aa_ce_s, axis_distance_aa_ce_e,
1460            axis_distance_aa_s_ce, axis_distance_aa_s_s, axis_distance_aa_s_e,
1461            axis_distance_aa_e_ce, axis_distance_aa_e_s, axis_distance_aa_e_e,
1462        );
1463        out
1464    }
1465
1466    /// Build a map from CID to human-readable constraint description.
1467    /// Must be called AFTER calc_jacobian/solve/compute_dof, which populates
1468    /// the `cid` field on each constraint instance.
1469    pub fn constraint_labels(&self) -> std::collections::HashMap<u32, String> {
1470        let mut m = std::collections::HashMap::new();
1471        for r in self.points.refs() { let p = &self.points[r]; m.insert(p.cid, format!("point:{}", p.name)); }
1472        for r in self.lines.refs() { let l = &self.lines[r]; m.insert(l.cid, format!("line:{}", l.name)); }
1473        for r in self.arcs.refs() { let a = &self.arcs[r]; m.insert(a.cid, format!("arc:{}", a.name)); }
1474        // Helper to get entity names
1475        let pn = |r: Ref<Point>| self.points.get(r).map(|p| p.name.as_str()).unwrap_or("?").to_string();
1476        let ln = |r: Ref<Line>| self.lines.get(r).map(|l| l.name.as_str()).unwrap_or("?").to_string();
1477        let an = |r: Ref<Arc>| self.arcs.get(r).map(|a| a.name.as_str()).unwrap_or("?").to_string();
1478        for c in &self.coincident_pp { m.insert(c.cid, format!("coinc_pp:{},{}", pn(c.a), pn(c.b))); }
1479        for c in &self.coincident_lp1 { m.insert(c.cid, format!("coinc_lp1:{},{}", ln(c.line), pn(c.point))); }
1480        for c in &self.coincident_lp2 { m.insert(c.cid, format!("coinc_lp2:{},{}", ln(c.line), pn(c.point))); }
1481        for c in &self.coincident_ll11 { m.insert(c.cid, format!("coinc_ll11:{},{}", ln(c.a), ln(c.b))); }
1482        for c in &self.coincident_ll12 { m.insert(c.cid, format!("coinc_ll12:{},{}", ln(c.a), ln(c.b))); }
1483        for c in &self.coincident_ll21 { m.insert(c.cid, format!("coinc_ll21:{},{}", ln(c.a), ln(c.b))); }
1484        for c in &self.coincident_ll22 { m.insert(c.cid, format!("coinc_ll22:{},{}", ln(c.a), ln(c.b))); }
1485        for c in &self.distance_pp { m.insert(c.cid, format!("dist_pp:{},{}", pn(c.a), pn(c.b))); }
1486        for c in &self.hdistance_pp { m.insert(c.cid, format!("hdist_pp:{},{}", pn(c.a), pn(c.b))); }
1487        for c in &self.vdistance_pp { m.insert(c.cid, format!("vdist_pp:{},{}", pn(c.a), pn(c.b))); }
1488        for c in &self.point_on_line { m.insert(c.cid, format!("on_line:{},{}", pn(c.point), ln(c.line))); }
1489        for c in &self.midpoint { m.insert(c.cid, format!("midpoint:{},{}", pn(c.point), ln(c.line))); }
1490        for c in &self.point_on_arc { m.insert(c.cid, format!("on_arc:{},{}", pn(c.point), an(c.arc))); }
1491        for c in &self.parallel { m.insert(c.cid, format!("parallel:{},{}", ln(c.a), ln(c.b))); }
1492        for c in &self.arc_line_parallel { m.insert(c.cid, format!("parallel:{},{}", an(c.arc), ln(c.line))); }
1493        for c in &self.arc_arc_parallel { m.insert(c.cid, format!("parallel:{},{}", an(c.a), an(c.b))); }
1494        for c in &self.perpendicular { m.insert(c.cid, format!("perp:{},{}", ln(c.a), ln(c.b))); }
1495        for c in &self.collinear { m.insert(c.cid, format!("collinear:{},{}", ln(c.a), ln(c.b))); }
1496        for c in &self.equal_length { m.insert(c.cid, format!("eq_len:{},{}", ln(c.a), ln(c.b))); }
1497        for c in &self.angle { m.insert(c.cid, format!("angle:{},{}", ln(c.a), ln(c.b))); }
1498        for c in &self.tangent_la { m.insert(c.cid, format!("tangent_la:{},{}", ln(c.line), an(c.arc))); }
1499        for c in &self.coincident_lp2_arc_start { m.insert(c.cid, format!("coinc_lp2_as:{},{}", ln(c.line), an(c.arc))); }
1500        for c in &self.tangent_aa { m.insert(c.cid, format!("tangent_aa:{},{}", an(c.a), an(c.b))); }
1501        for c in &self.concentric { m.insert(c.cid, format!("concentric:{},{}", an(c.a), an(c.b))); }
1502        for c in &self.equal_radius { m.insert(c.cid, format!("eq_radius:{},{}", an(c.a), an(c.b))); }
1503        for c in &self.coincident_arc_center { m.insert(c.cid, format!("coinc_arc_center:{},{}", pn(c.point), an(c.arc))); }
1504        for c in &self.coincident_arc_start { m.insert(c.cid, format!("coinc_arc_start:{},{}", pn(c.point), an(c.arc))); }
1505        for c in &self.coincident_arc_end { m.insert(c.cid, format!("coinc_arc_end:{},{}", pn(c.point), an(c.arc))); }
1506        for c in &self.coincident_lp1_arc_center { m.insert(c.cid, format!("coinc_lp1_ac:{},{}", ln(c.line), an(c.arc))); }
1507        for c in &self.coincident_lp2_arc_center { m.insert(c.cid, format!("coinc_lp2_ac:{},{}", ln(c.line), an(c.arc))); }
1508        for c in &self.coincident_lp1_arc_start { m.insert(c.cid, format!("coinc_lp1_as:{},{}", ln(c.line), an(c.arc))); }
1509        for c in &self.coincident_lp2_arc_start { m.insert(c.cid, format!("coinc_lp2_as:{},{}", ln(c.line), an(c.arc))); }
1510        for c in &self.coincident_lp1_arc_end { m.insert(c.cid, format!("coinc_lp1_ae:{},{}", ln(c.line), an(c.arc))); }
1511        for c in &self.coincident_lp2_arc_end { m.insert(c.cid, format!("coinc_lp2_ae:{},{}", ln(c.line), an(c.arc))); }
1512        for c in &self.coincident_arc_center_start { m.insert(c.cid, format!("coinc_ac_as:{},{}", an(c.a), an(c.b))); }
1513        for c in &self.coincident_arc_center_end { m.insert(c.cid, format!("coinc_ac_ae:{},{}", an(c.a), an(c.b))); }
1514        for c in &self.coincident_arc_start_center { m.insert(c.cid, format!("coinc_as_ac:{},{}", an(c.a), an(c.b))); }
1515        for c in &self.coincident_arc_end_center { m.insert(c.cid, format!("coinc_ae_ac:{},{}", an(c.a), an(c.b))); }
1516        for c in &self.coincident_arc_start_start { m.insert(c.cid, format!("coinc_as_as:{},{}", an(c.a), an(c.b))); }
1517        for c in &self.coincident_arc_start_end { m.insert(c.cid, format!("coinc_as_ae:{},{}", an(c.a), an(c.b))); }
1518        for c in &self.coincident_arc_end_start { m.insert(c.cid, format!("coinc_ae_as:{},{}", an(c.a), an(c.b))); }
1519        for c in &self.coincident_arc_end_end { m.insert(c.cid, format!("coinc_ae_ae:{},{}", an(c.a), an(c.b))); }
1520        for c in &self.line_p1_on_line { m.insert(c.cid, format!("lp1_on_l:{},{}", ln(c.a), ln(c.b))); }
1521        for c in &self.line_p2_on_line { m.insert(c.cid, format!("lp2_on_l:{},{}", ln(c.a), ln(c.b))); }
1522        for c in &self.line_p1_on_arc { m.insert(c.cid, format!("lp1_on_a:{},{}", ln(c.line), an(c.arc))); }
1523        for c in &self.line_p2_on_arc { m.insert(c.cid, format!("lp2_on_a:{},{}", ln(c.line), an(c.arc))); }
1524        for c in &self.distance_pl { m.insert(c.cid, format!("dist_pl:{},{}", pn(c.point), ln(c.line))); }
1525        for c in &self.distance_lp1l { m.insert(c.cid, format!("dist_lp1l:{},{}", ln(c.a), ln(c.b))); }
1526        for c in &self.distance_lp2l { m.insert(c.cid, format!("dist_lp2l:{},{}", ln(c.a), ln(c.b))); }
1527        for c in &self.distance_arc_center_l { m.insert(c.cid, format!("dist_ac_l:{},{}", an(c.arc), ln(c.line))); }
1528        for c in &self.distance_arc_start_l { m.insert(c.cid, format!("dist_as_l:{},{}", an(c.arc), ln(c.line))); }
1529        for c in &self.distance_arc_end_l { m.insert(c.cid, format!("dist_ae_l:{},{}", an(c.arc), ln(c.line))); }
1530        for c in &self.distance_concentric { m.insert(c.cid, format!("dist_concentric:{},{}", an(c.a), an(c.b))); }
1531        for c in &self.symmetry_ll { m.insert(c.cid, format!("sym_ll:{},{},{}", ln(c.a), ln(c.b), ln(c.c))); }
1532        for c in &self.symmetry_pp { m.insert(c.cid, format!("sym_pp:{},{},{}", pn(c.a), ln(c.line), pn(c.c))); }
1533        for c in &self.symmetry_aa { m.insert(c.cid, format!("sym_aa:{},{},{}", an(c.a), ln(c.line), an(c.c))); }
1534        for ec in &self.expr_constraints { m.insert(ec.cid, format!("expr:{}", ec.description)); }
1535        for (i, ec) in self.expr_constraints.iter().enumerate() {
1536            // expr_constraints are assigned CIDs starting after all macro-generated ones
1537            // but we don't have easy access to that base — approximate with description
1538            let _ = i;
1539            let label = format!("expr:{}", ec.description);
1540            // We don't know the CID here — it would need to be tracked. Skip for now.
1541            let _ = label;
1542        }
1543        m
1544    }
1545
1546    /// Fix up arc fields after loading old files that lack radius_b/rotation.
1547    /// Detects the sentinel default (radius_b == 0.0) and sets radius_b to
1548    /// match radius, rotation to 0, as optimizable params.
1549    pub fn fixup_after_load(&mut self) {
1550        let refs: Vec<_> = self.arcs.refs().collect();
1551        for r in refs {
1552            if self.arcs[r].radius_b.value == 0.0 && !self.arcs[r].is_ellipse {
1553                let rv = self.arcs[r].radius.value;
1554                self.arcs[r].radius_b = Param::new(rv);
1555                self.arcs[r].rotation = Param::fixed(0.0);
1556            }
1557        }
1558    }
1559
1560    /// Add an expression constraint. The expression should evaluate to 0
1561    /// when satisfied. Symbol resolution and differentiation happen at
1562    /// solve() time.
1563    pub fn add_expr_constraint(&mut self, expr: arael_sym::E, description: String) {
1564        self.expr_constraints.push(ExpressionConstraint::new_unresolved(expr, description));
1565    }
1566
1567    /// Rebuild expr_constraints from dimensions that have expr_str.
1568    /// Called at the start of every solve() since the set of optimizable
1569    /// params can change between solves (lock/unlock).
1570    fn rebuild_expr_constraints(&mut self) {
1571        self.expr_constraints.clear();
1572        let has_expr = self.dimensions.iter().any(|d| d.expr_str.is_some());
1573        let has_user_params = !self.user_params.is_empty();
1574        let has_range = self.dimensions.iter().any(|d| d.range.is_some());
1575        if !has_expr && !has_user_params && !has_range {
1576            for d in &mut self.dimensions { d.broken = false; }
1577            return;
1578        }
1579
1580        // Reset broken flags so SymbolBag always starts fresh --
1581        // stale flags from a previous solve can hide circular refs.
1582        for d in &mut self.dimensions { d.broken = false; }
1583        for p in &mut self.user_params { p.broken = false; }
1584
1585        // Need param indices assigned for SymbolBag; serialize to assign them.
1586        {
1587            let mut tmp = std::vec::Vec::new();
1588            self.serialize64(&mut tmp);
1589        }
1590        let mut bag = SymbolBag::build(self);
1591
1592        // Detect broken user params first (they feed into dimensions).
1593        // Process in order so earlier params that break get frozen before
1594        // downstream params/dims are checked.
1595        for i in 0..self.user_params.len() {
1596            let expr_str = &self.user_params[i].expr_str;
1597            // Pure numeric literals are never broken
1598            if expr_str.trim().parse::<f64>().is_ok() { continue; }
1599            let is_broken = if let Ok(parsed) = arael_sym::parse(expr_str) {
1600                let expanded = expr_constraint::expand_derived(&parsed, &bag);
1601                !expanded.symbols().iter().all(|sym|
1602                    bag.param_indices.contains_key(sym.as_str())
1603                    || bag.dim_values.contains_key(sym.as_str())
1604                )
1605            } else {
1606                true
1607            };
1608            self.user_params[i].broken = is_broken;
1609            if is_broken {
1610                bag.derived.remove(&self.user_params[i].name);
1611                bag.dim_values.insert(
1612                    self.user_params[i].name.clone(),
1613                    self.user_params[i].value,
1614                );
1615            }
1616        }
1617
1618        // Detect broken references and create expression constraints.
1619        // Process in order so broken dims get frozen in the bag before
1620        // downstream dims that reference them are checked.
1621        for i in 0..self.dimensions.len() {
1622            // Derived dimensions don't create constraints
1623            if self.dimensions[i].derived { continue; }
1624            if let Some(ref expr_str) = self.dimensions[i].expr_str {
1625                let is_broken = if let Ok(parsed) = arael_sym::parse(expr_str) {
1626                    let expanded = expr_constraint::expand_derived(&parsed, &bag);
1627                    let all_resolved = expanded.symbols().iter().all(|sym|
1628                        bag.param_indices.contains_key(sym.as_str())
1629                        || bag.dim_values.contains_key(sym.as_str())
1630                    );
1631                    if all_resolved {
1632                        // Normal: measured - expr = 0
1633                        let measured = self.dimensions[i].measured_symbol(self);
1634                        let residual = measured - parsed;
1635                        let desc = format!("{} = {}", self.dimensions[i].name, expr_str);
1636                        self.expr_constraints.push(
1637                            ExpressionConstraint::new_unresolved(residual, desc));
1638                        false
1639                    } else {
1640                        true
1641                    }
1642                } else {
1643                    true
1644                };
1645
1646                self.dimensions[i].broken = is_broken;
1647                if is_broken {
1648                    // Freeze in bag so downstream dims see a constant
1649                    bag.derived.remove(&self.dimensions[i].name);
1650                    bag.dim_values.insert(
1651                        self.dimensions[i].name.clone(),
1652                        self.dimensions[i].value,
1653                    );
1654                    // Fallback: constrain to last computed value
1655                    let measured = self.dimensions[i].measured_symbol(self);
1656                    let residual = measured - arael_sym::constant(self.dimensions[i].value);
1657                    let desc = format!("{} = {} [broken]", self.dimensions[i].name, self.dimensions[i].value);
1658                    self.expr_constraints.push(
1659                        ExpressionConstraint::new_unresolved(residual, desc));
1660                }
1661            } else {
1662                self.dimensions[i].broken = false;
1663            }
1664        }
1665
1666        // Range dimensions: barrier residuals (piecewise-zero inside the
1667        // feasible region). Each contributes one direct ExpressionConstraint
1668        // whose `expr` IS the residual (not measured - expr, as in the
1669        // equality path above). Live bounds re-parse + resolve against the
1670        // current SymbolBag; a bound with unresolved free symbols marks
1671        // the dimension broken (same cascade as expression dims above).
1672        for i in 0..self.dimensions.len() {
1673            let rb = match self.dimensions[i].range.clone() {
1674                Some(rb) => rb,
1675                None => continue,
1676            };
1677            let resolve_value = |rv: &dimensions::RangeValue, bag: &SymbolBag|
1678                -> Option<arael_sym::E>
1679            {
1680                match rv {
1681                    dimensions::RangeValue::Literal(v) => Some(arael_sym::constant(*v)),
1682                    dimensions::RangeValue::Live(src) => {
1683                        let parsed = arael_sym::parse(src).ok()?;
1684                        let expanded = expr_constraint::expand_derived(&parsed, bag);
1685                        let all_resolved = expanded.symbols().iter().all(|sym|
1686                            bag.param_indices.contains_key(sym.as_str())
1687                            || bag.dim_values.contains_key(sym.as_str()));
1688                        if all_resolved { Some(expanded) } else { None }
1689                    }
1690                }
1691            };
1692            let resolved = match &rb {
1693                dimensions::RangeBound::Min(v) =>
1694                    resolve_value(v, &bag).map(dimensions::ResolvedBound::Min),
1695                dimensions::RangeBound::Max(v) =>
1696                    resolve_value(v, &bag).map(dimensions::ResolvedBound::Max),
1697                dimensions::RangeBound::Between(lo, hi) => {
1698                    match (resolve_value(lo, &bag), resolve_value(hi, &bag)) {
1699                        (Some(l), Some(h)) => Some(dimensions::ResolvedBound::Between(l, h)),
1700                        _ => None,
1701                    }
1702                }
1703            };
1704            let Some(resolved) = resolved else {
1705                self.dimensions[i].broken = true;
1706                continue;
1707            };
1708            self.dimensions[i].broken = false;
1709            let measured = self.dimensions[i].measured_symbol(self);
1710            let residual = Dimension::range_residual(&resolved, measured);
1711            let desc = format!("{} {}", self.dimensions[i].name,
1712                match &rb {
1713                    dimensions::RangeBound::Min(v) => format!(">= {}", v),
1714                    dimensions::RangeBound::Max(v) => format!("<= {}", v),
1715                    dimensions::RangeBound::Between(lo, hi) => format!("in {} to {}", lo, hi),
1716                });
1717            let mut ec = ExpressionConstraint::new_unresolved(residual, desc);
1718            // Range dimensions are one-sided barriers: the residual is
1719            // zero inside the feasible band and linear outside it. Tag
1720            // the label so DOF rank detection can strip them -- their
1721            // Jacobian row would otherwise flip in/out of rank as the
1722            // geometry crosses the bound, reporting DOF 2 at the bound
1723            // and DOF 3 off it for the same sketch.
1724            ec.label = "range";
1725            self.expr_constraints.push(ec);
1726        }
1727    }
1728
1729    /// Validate an expression string: parse it and check all symbols resolve.
1730    /// Returns Err with a description if invalid.
1731    pub fn validate_expr(&mut self, expr_str: &str) -> Result<(), String> {
1732        let parsed = arael_sym::parse(expr_str).map_err(|e| e.to_string())?;
1733        {
1734            let mut tmp = std::vec::Vec::new();
1735            self.serialize64(&mut tmp);
1736        }
1737        let bag = SymbolBag::build(self);
1738        let expanded = expr_constraint::expand_derived(&parsed, &bag);
1739        let unresolved: Vec<String> = expanded.symbols().into_iter().filter(|sym|
1740            !bag.param_indices.contains_key(sym.as_str())
1741            && !bag.dim_values.contains_key(sym.as_str())
1742        ).collect();
1743        if !unresolved.is_empty() {
1744            return Err(format!("Unknown symbol: {}", unresolved.join(", ")));
1745        }
1746        Ok(())
1747    }
1748
1749    /// Validate a user parameter name. Returns Err if the name is empty,
1750    /// a duplicate, a system name pattern, or already used by an entity.
1751    pub fn validate_param_name(&self, name: &str, exclude_index: Option<usize>) -> Result<(), String> {
1752        let name = name.trim();
1753        if name.is_empty() {
1754            return Err("Name cannot be empty".into());
1755        }
1756        // Must be a valid identifier: alphanumeric + underscore, not starting with digit
1757        if name.bytes().next().is_none_or(|b| b.is_ascii_digit()) {
1758            return Err("Name cannot start with a digit".into());
1759        }
1760        if !name.bytes().all(|b| b.is_ascii_alphanumeric() || b == b'_') {
1761            return Err("Name can only contain letters, digits, and underscores".into());
1762        }
1763        // Not a system name pattern (d0, L0, P0, A0, etc.)
1764        if is_system_name(name) {
1765            return Err(format!("'{}' conflicts with system naming (d/L/P/A + number)", name));
1766        }
1767        // Not a duplicate of another user param
1768        for (i, p) in self.user_params.iter().enumerate() {
1769            if Some(i) == exclude_index { continue; }
1770            if p.name == name {
1771                return Err(format!("Parameter '{}' already exists", name));
1772            }
1773        }
1774        // Not a dimension name
1775        for d in &self.dimensions {
1776            if d.name == name {
1777                return Err(format!("'{}' is already used by a dimension", name));
1778            }
1779        }
1780        Ok(())
1781    }
1782
1783    /// Get a display-friendly name for a point. For helper points, resolves
1784    /// through bridge constraints to show the semantic origin (e.g. "L0.p1").
1785    pub fn point_display_name(&self, r: Ref<Point>) -> String {
1786        let p = &self.points[r];
1787        if !p.helper { return p.name.clone(); }
1788        for c in &self.coincident_lp1 {
1789            if c.point == r { return format!("{}.p1", self.lines[c.line].name); }
1790        }
1791        for c in &self.coincident_lp2 {
1792            if c.point == r { return format!("{}.p2", self.lines[c.line].name); }
1793        }
1794        for c in &self.coincident_arc_center {
1795            if c.point == r { return format!("{}.center", self.arcs[c.arc].name); }
1796        }
1797        for c in &self.coincident_arc_start {
1798            if c.point == r { return format!("{}.start", self.arcs[c.arc].name); }
1799        }
1800        for c in &self.coincident_arc_end {
1801            if c.point == r { return format!("{}.end", self.arcs[c.arc].name); }
1802        }
1803        p.name.clone()
1804    }
1805
1806    /// Add an expression-based dimension. The expression string is parsed
1807    /// Build a `cid -> dimension_name` map for constraints that are
1808    /// backed by a dimension (`d<n>`). Deleting the dimension via
1809    /// `delete d<n>` removes both the dimension and its backing
1810    /// constraint, which is the user-facing way to edit or remove
1811    /// these constraints -- `C<nid>` is visible in `list` but is not
1812    /// a valid handle for `delete` (see [`tools::find_constraint_by_name`]).
1813    ///
1814    /// Coverage focuses on the dimension kinds that write into
1815    /// constraint collections: `HDistance`, `VDistance`,
1816    /// `PointPointDistance`, `PointLineDistance`, `Angle`,
1817    /// `ConcentricDistance`, `LineLineDistance`. Dimensions that only
1818    /// feed expression constraints (`LineLength`, `ArcRadius`,
1819    /// `ArcRadiusB`, `ArcSweep`, `LineAngle`) don't have a backing
1820    /// cid in any collection and are omitted.
1821    pub fn dimension_cid_name_map(&self) -> std::collections::HashMap<u32, String> {
1822        use crate::dimensions::{DimensionEndpoint as E, DimensionKind as K};
1823        let mut out = std::collections::HashMap::new();
1824        for dim in &self.dimensions {
1825            // Each arm scans the collection(s) that the corresponding
1826            // add path pushes into. See `push_axis_distance` and
1827            // `add_dimension` in arael-sketch/src/actions.rs for the
1828            // inverse mapping.
1829            let cid = match &dim.kind {
1830                K::HDistance(a, b) => find_axis_cid(self, a, b, true),
1831                K::VDistance(a, b) => find_axis_cid(self, a, b, false),
1832                K::PointPointDistance(a, b) => find_distance_cid(self, a, b),
1833                K::PointLineDistance(ep, l) => find_point_line_cid(self, ep, *l),
1834                K::Angle(la, lb, _) => self.angle.iter()
1835                    .find(|c| c.a == *la && c.b == *lb)
1836                    .map(|c| c.cid),
1837                K::ConcentricDistance(a, b) => self.distance_concentric.iter()
1838                    .find(|c| (c.a == *a && c.b == *b) || (c.a == *b && c.b == *a))
1839                    .map(|c| c.cid),
1840                K::LineLineDistance(a, b) => {
1841                    // LineLineDistance uses point-to-line with LineP1(b) as anchor.
1842                    find_point_line_cid(self, &E::LineP1(*b), *a)
1843                }
1844                K::LineLength(_) | K::ArcRadius(_) | K::ArcRadiusB(_)
1845                    | K::ArcSweep(_) | K::LineAngle(_) | K::ArcRotation(_) => None,
1846            };
1847            if let Some(cid) = cid {
1848                out.insert(cid, dim.name.clone());
1849            }
1850        }
1851        // Expression-backed dimensions (LineLength, ArcRadius,
1852        // ArcRadiusB, ArcSweep, LineAngle, plus any dimension with an
1853        // expr_str or range bound) feed into self.expr_constraints.
1854        // Their descriptions always start with the dimension name
1855        // followed by a space -- parse that prefix to map cid -> name.
1856        for ec in &self.expr_constraints {
1857            if let Some((dim_name, _)) = ec.description.split_once(' ') {
1858                let is_dim = dim_name.starts_with('d')
1859                    && dim_name.len() > 1
1860                    && dim_name[1..].chars().all(|c| c.is_ascii_digit());
1861                if is_dim {
1862                    out.insert(ec.cid, dim_name.to_string());
1863                }
1864            }
1865        }
1866        out
1867    }
1868
1869    /// Look up a constraint by its user-visible name (e.g. `"C1"`,
1870    /// `"CL0H"`) and return the descriptive part of its
1871    /// `list_constraints` entry. Returns `None` if no active
1872    /// constraint bears that name.
1873    ///
1874    /// The search walks `list_constraints` and strips the `"<name>: "`
1875    /// prefix, so every constraint that appears in the canonical list
1876    /// output is findable -- including dimension-managed distance
1877    /// constraints that have a `C<nid>` name but no `ConstraintId`
1878    /// variant.
1879    pub fn find_constraint_description(&self, name: &str) -> Option<String> {
1880        let prefix = format!("{}: ", name);
1881        for line in self.list_constraints() {
1882            if let Some(rest) = line.strip_prefix(&prefix) {
1883                return Some(rest.to_string());
1884            }
1885        }
1886        None
1887    }
1888
1889    /// List all active constraints as human-readable strings.
1890    pub fn list_constraints(&self) -> Vec<String> {
1891        let mut out = Vec::new();
1892        // Entity-level flags
1893        for r in self.lines.refs() {
1894            let l = &self.lines[r];
1895            if l.constraints.horizontal { out.push(format!("{}: horizontal {}", format_flag_name(&l.name, 'H'), l.name)); }
1896            if l.constraints.vertical { out.push(format!("{}: vertical {}", format_flag_name(&l.name, 'V'), l.name)); }
1897            if l.constraints.has_length { out.push(format!("length {} = {}", l.name, l.constraints.length)); }
1898            if l.constraints.has_angle { out.push(format!("xangle {} = {:.4}", l.name, l.constraints.target_angle.to_degrees())); }
1899            if !l.p1.optimize { out.push(format!("lock {}.p1", l.name)); }
1900            if !l.p2.optimize { out.push(format!("lock {}.p2", l.name)); }
1901        }
1902        for r in self.points.refs() {
1903            let p = &self.points[r];
1904            if p.constraints.has_fix_x || p.constraints.has_fix_y {
1905                out.push(format!("lock {}", p.name));
1906            }
1907        }
1908        for r in self.arcs.refs() {
1909            let a = &self.arcs[r];
1910            if a.constraints.has_target_radius { out.push(format!("radius {} = {}", a.name, a.constraints.target_radius)); }
1911            if a.constraints.has_target_sweep { out.push(format!("sweep {} = {:.2} deg", a.name, a.constraints.target_sweep.to_degrees())); }
1912            if !a.center.optimize { out.push(format!("lock {}.center", a.name)); }
1913        }
1914        // Constraint vectors (use macros to reduce repetition)
1915        macro_rules! list_cross {
1916            ($vec:expr, $name:literal, $fa:ident, $fb:ident) => {
1917                for c in &$vec {
1918                    let na = &self.lines[c.$fa].name;
1919                    let nb = &self.lines[c.$fb].name;
1920                    out.push(format!("C{}: {} {} {}", c.nid, $name, na, nb));
1921                }
1922            };
1923        }
1924        list_cross!(self.parallel, "parallel", a, b);
1925        for c in &self.arc_line_parallel {
1926            let na = &self.arcs[c.arc].name;
1927            let nb = &self.lines[c.line].name;
1928            out.push(format!("C{}: parallel {} {}", c.nid, na, nb));
1929        }
1930        for c in &self.arc_arc_parallel {
1931            let na = &self.arcs[c.a].name;
1932            let nb = &self.arcs[c.b].name;
1933            out.push(format!("C{}: parallel {} {}", c.nid, na, nb));
1934        }
1935        list_cross!(self.perpendicular, "perpendicular", a, b);
1936        list_cross!(self.collinear, "collinear", a, b);
1937        list_cross!(self.equal_length, "equal", a, b);
1938        // Coincident: suppress bridge constraints (helper point bridges)
1939        for c in &self.coincident_pp {
1940            if !self.points[c.a].helper && !self.points[c.b].helper {
1941                out.push(format!("C{}: coincident {} {}", c.nid, self.points[c.a].name, self.points[c.b].name));
1942            }
1943        }
1944        for c in &self.coincident_ll11 { out.push(format!("C{}: coincident {}.p1 {}.p1", c.nid, self.lines[c.a].name, self.lines[c.b].name)); }
1945        for c in &self.coincident_ll12 { out.push(format!("C{}: coincident {}.p1 {}.p2", c.nid, self.lines[c.a].name, self.lines[c.b].name)); }
1946        for c in &self.coincident_ll21 { out.push(format!("C{}: coincident {}.p2 {}.p1", c.nid, self.lines[c.a].name, self.lines[c.b].name)); }
1947        for c in &self.coincident_ll22 { out.push(format!("C{}: coincident {}.p2 {}.p2", c.nid, self.lines[c.a].name, self.lines[c.b].name)); }
1948        for c in &self.coincident_lp1 {
1949            if !self.points[c.point].helper {
1950                out.push(format!("C{}: coincident {}.p1 {}", c.nid, self.lines[c.line].name, self.points[c.point].name));
1951            }
1952        }
1953        for c in &self.coincident_lp2 {
1954            if !self.points[c.point].helper {
1955                out.push(format!("C{}: coincident {}.p2 {}", c.nid, self.lines[c.line].name, self.points[c.point].name));
1956            }
1957        }
1958        // Point-Arc coincident: suppress helper bridges
1959        for c in &self.coincident_arc_center {
1960            if !self.points[c.point].helper {
1961                out.push(format!("C{}: coincident {} {}.center", c.nid, self.points[c.point].name, self.arcs[c.arc].name));
1962            }
1963        }
1964        for c in &self.coincident_arc_start {
1965            if !self.points[c.point].helper {
1966                out.push(format!("C{}: coincident {} {}.start", c.nid, self.points[c.point].name, self.arcs[c.arc].name));
1967            }
1968        }
1969        for c in &self.coincident_arc_end {
1970            if !self.points[c.point].helper {
1971                out.push(format!("C{}: coincident {} {}.end", c.nid, self.points[c.point].name, self.arcs[c.arc].name));
1972            }
1973        }
1974        // Line-Arc coincident
1975        for c in &self.coincident_lp1_arc_center { out.push(format!("C{}: coincident {}.p1 {}.center", c.nid, self.lines[c.line].name, self.arcs[c.arc].name)); }
1976        for c in &self.coincident_lp2_arc_center { out.push(format!("C{}: coincident {}.p2 {}.center", c.nid, self.lines[c.line].name, self.arcs[c.arc].name)); }
1977        for c in &self.coincident_lp1_arc_start { out.push(format!("C{}: coincident {}.p1 {}.start", c.nid, self.lines[c.line].name, self.arcs[c.arc].name)); }
1978        for c in &self.coincident_lp2_arc_start { out.push(format!("C{}: coincident {}.p2 {}.start", c.nid, self.lines[c.line].name, self.arcs[c.arc].name)); }
1979        for c in &self.coincident_lp1_arc_end { out.push(format!("C{}: coincident {}.p1 {}.end", c.nid, self.lines[c.line].name, self.arcs[c.arc].name)); }
1980        for c in &self.coincident_lp2_arc_end { out.push(format!("C{}: coincident {}.p2 {}.end", c.nid, self.lines[c.line].name, self.arcs[c.arc].name)); }
1981        // Arc-Arc coincident
1982        for c in &self.coincident_arc_center_start { out.push(format!("C{}: coincident {}.center {}.start", c.nid, self.arcs[c.a].name, self.arcs[c.b].name)); }
1983        for c in &self.coincident_arc_center_end { out.push(format!("C{}: coincident {}.center {}.end", c.nid, self.arcs[c.a].name, self.arcs[c.b].name)); }
1984        for c in &self.coincident_arc_start_center { out.push(format!("C{}: coincident {}.start {}.center", c.nid, self.arcs[c.a].name, self.arcs[c.b].name)); }
1985        for c in &self.coincident_arc_end_center { out.push(format!("C{}: coincident {}.end {}.center", c.nid, self.arcs[c.a].name, self.arcs[c.b].name)); }
1986        for c in &self.coincident_arc_start_start { out.push(format!("C{}: coincident {}.start {}.start", c.nid, self.arcs[c.a].name, self.arcs[c.b].name)); }
1987        for c in &self.coincident_arc_start_end { out.push(format!("C{}: coincident {}.start {}.end", c.nid, self.arcs[c.a].name, self.arcs[c.b].name)); }
1988        for c in &self.coincident_arc_end_start { out.push(format!("C{}: coincident {}.end {}.start", c.nid, self.arcs[c.a].name, self.arcs[c.b].name)); }
1989        for c in &self.coincident_arc_end_end { out.push(format!("C{}: coincident {}.end {}.end", c.nid, self.arcs[c.a].name, self.arcs[c.b].name)); }
1990        // Line endpoint on line/arc
1991        for c in &self.line_p1_on_line { out.push(format!("C{}: point_on {}.p1 {}", c.nid, self.lines[c.a].name, self.lines[c.b].name)); }
1992        for c in &self.line_p2_on_line { out.push(format!("C{}: point_on {}.p2 {}", c.nid, self.lines[c.a].name, self.lines[c.b].name)); }
1993        for c in &self.line_p1_on_arc { out.push(format!("C{}: point_on {}.p1 {}", c.nid, self.lines[c.line].name, self.arcs[c.arc].name)); }
1994        for c in &self.line_p2_on_arc { out.push(format!("C{}: point_on {}.p2 {}", c.nid, self.lines[c.line].name, self.arcs[c.arc].name)); }
1995        for c in &self.angle { out.push(format!("C{}: angle {} {} = {:.1}deg", c.nid, self.lines[c.a].name, self.lines[c.b].name, c.angle.to_degrees())); }
1996        for c in &self.tangent_la { out.push(format!("C{}: tangent {} {}", c.nid, self.lines[c.line].name, self.arcs[c.arc].name)); }
1997        for c in &self.tangent_aa { out.push(format!("C{}: tangent {} {}", c.nid, self.arcs[c.a].name, self.arcs[c.b].name)); }
1998        for c in &self.concentric { out.push(format!("C{}: concentric {} {}", c.nid, self.arcs[c.a].name, self.arcs[c.b].name)); }
1999        for c in &self.equal_radius { out.push(format!("C{}: equal {} {}", c.nid, self.arcs[c.a].name, self.arcs[c.b].name)); }
2000        for c in &self.symmetry_ll { out.push(format!("C{}: symmetry {} {} {}", c.nid, self.lines[c.a].name, self.lines[c.b].name, self.lines[c.c].name)); }
2001        for c in &self.symmetry_pp { out.push(format!("C{}: symmetry {} {} {}", c.nid, self.point_display_name(c.a), self.lines[c.line].name, self.point_display_name(c.c))); }
2002        for c in &self.symmetry_aa { out.push(format!("C{}: symmetry {} {} {}", c.nid, self.arcs[c.a].name, self.lines[c.line].name, self.arcs[c.c].name)); }
2003        // Point-based constraints: use display names to resolve helpers
2004        for c in &self.point_on_line { out.push(format!("C{}: point_on {} {}", c.nid, self.point_display_name(c.point), self.lines[c.line].name)); }
2005        for c in &self.point_on_arc { out.push(format!("C{}: point_on {} {}", c.nid, self.point_display_name(c.point), self.arcs[c.arc].name)); }
2006        for c in &self.midpoint { out.push(format!("C{}: midpoint {} {}", c.nid, self.point_display_name(c.point), self.lines[c.line].name)); }
2007        for c in &self.midpoint_lp1 { out.push(format!("C{}: midpoint {}.p1 {}", c.nid, self.lines[c.line].name, self.lines[c.target].name)); }
2008        for c in &self.midpoint_lp2 { out.push(format!("C{}: midpoint {}.p2 {}", c.nid, self.lines[c.line].name, self.lines[c.target].name)); }
2009        for c in &self.midpoint_arc_start { out.push(format!("C{}: midpoint {}.start {}", c.nid, self.arcs[c.arc].name, self.lines[c.line].name)); }
2010        for c in &self.midpoint_arc_end { out.push(format!("C{}: midpoint {}.end {}", c.nid, self.arcs[c.arc].name, self.lines[c.line].name)); }
2011        for c in &self.midpoint_arc_point { out.push(format!("C{}: midpoint {} {}", c.nid, self.point_display_name(c.point), self.arcs[c.arc].name)); }
2012        for c in &self.midpoint_lp1_arc { out.push(format!("C{}: midpoint {}.p1 {}", c.nid, self.lines[c.line].name, self.arcs[c.arc].name)); }
2013        for c in &self.midpoint_lp2_arc { out.push(format!("C{}: midpoint {}.p2 {}", c.nid, self.lines[c.line].name, self.arcs[c.arc].name)); }
2014        for c in &self.midpoint_arc_start_arc { out.push(format!("C{}: midpoint {}.start {}", c.nid, self.arcs[c.a].name, self.arcs[c.b].name)); }
2015        for c in &self.midpoint_arc_end_arc { out.push(format!("C{}: midpoint {}.end {}", c.nid, self.arcs[c.a].name, self.arcs[c.b].name)); }
2016        for c in &self.distance_pp { out.push(format!("C{}: distance {} {} = {}", c.nid, self.point_display_name(c.a), self.point_display_name(c.b), c.distance.abs())); }
2017        for c in &self.hdistance_pp { out.push(format!("C{}: hdistance {} {} = {}", c.nid, self.point_display_name(c.a), self.point_display_name(c.b), c.distance.abs())); }
2018        for c in &self.vdistance_pp { out.push(format!("C{}: vdistance {} {} = {}", c.nid, self.point_display_name(c.a), self.point_display_name(c.b), c.distance.abs())); }
2019        for c in &self.distance_pl { out.push(format!("C{}: distance {} {} = {}", c.nid, self.point_display_name(c.point), self.lines[c.line].name, c.distance.abs())); }
2020        for c in &self.distance_lp1l { out.push(format!("C{}: distance {}.p1 {} = {}", c.nid, self.lines[c.a].name, self.lines[c.b].name, c.distance.abs())); }
2021        for c in &self.distance_lp2l { out.push(format!("C{}: distance {}.p2 {} = {}", c.nid, self.lines[c.a].name, self.lines[c.b].name, c.distance.abs())); }
2022        for c in &self.distance_arc_center_l { out.push(format!("C{}: distance {}.center {} = {}", c.nid, self.arcs[c.arc].name, self.lines[c.line].name, c.distance.abs())); }
2023        for c in &self.distance_arc_start_l { out.push(format!("C{}: distance {}.start {} = {}", c.nid, self.arcs[c.arc].name, self.lines[c.line].name, c.distance.abs())); }
2024        for c in &self.distance_arc_end_l { out.push(format!("C{}: distance {}.end {} = {}", c.nid, self.arcs[c.arc].name, self.lines[c.line].name, c.distance.abs())); }
2025        // Line-endpoint distance constraints
2026        for c in &self.distance_ll11 { out.push(format!("C{}: distance {}.p1 {}.p1 = {}", c.nid, self.lines[c.a].name, self.lines[c.b].name, c.distance.abs())); }
2027        for c in &self.distance_ll12 { out.push(format!("C{}: distance {}.p1 {}.p2 = {}", c.nid, self.lines[c.a].name, self.lines[c.b].name, c.distance.abs())); }
2028        for c in &self.distance_ll21 { out.push(format!("C{}: distance {}.p2 {}.p1 = {}", c.nid, self.lines[c.a].name, self.lines[c.b].name, c.distance.abs())); }
2029        for c in &self.distance_ll22 { out.push(format!("C{}: distance {}.p2 {}.p2 = {}", c.nid, self.lines[c.a].name, self.lines[c.b].name, c.distance.abs())); }
2030        for c in &self.distance_lp1 { out.push(format!("C{}: distance {}.p1 {} = {}", c.nid, self.lines[c.line].name, self.point_display_name(c.point), c.distance.abs())); }
2031        for c in &self.distance_lp2 { out.push(format!("C{}: distance {}.p2 {} = {}", c.nid, self.lines[c.line].name, self.point_display_name(c.point), c.distance.abs())); }
2032        for c in &self.distance_arc_center_p { out.push(format!("C{}: distance {}.center {} = {}", c.nid, self.arcs[c.arc].name, self.point_display_name(c.point), c.distance.abs())); }
2033        for c in &self.distance_arc_start_p { out.push(format!("C{}: distance {}.start {} = {}", c.nid, self.arcs[c.arc].name, self.point_display_name(c.point), c.distance.abs())); }
2034        for c in &self.distance_arc_end_p { out.push(format!("C{}: distance {}.end {} = {}", c.nid, self.arcs[c.arc].name, self.point_display_name(c.point), c.distance.abs())); }
2035        for c in &self.distance_arc_center_l1 { out.push(format!("C{}: distance {}.center {}.p1 = {}", c.nid, self.arcs[c.arc].name, self.lines[c.line].name, c.distance.abs())); }
2036        for c in &self.distance_arc_center_l2 { out.push(format!("C{}: distance {}.center {}.p2 = {}", c.nid, self.arcs[c.arc].name, self.lines[c.line].name, c.distance.abs())); }
2037        for c in &self.distance_arc_start_l1 { out.push(format!("C{}: distance {}.start {}.p1 = {}", c.nid, self.arcs[c.arc].name, self.lines[c.line].name, c.distance.abs())); }
2038        for c in &self.distance_arc_start_l2 { out.push(format!("C{}: distance {}.start {}.p2 = {}", c.nid, self.arcs[c.arc].name, self.lines[c.line].name, c.distance.abs())); }
2039        for c in &self.distance_arc_end_l1 { out.push(format!("C{}: distance {}.end {}.p1 = {}", c.nid, self.arcs[c.arc].name, self.lines[c.line].name, c.distance.abs())); }
2040        for c in &self.distance_arc_end_l2 { out.push(format!("C{}: distance {}.end {}.p2 = {}", c.nid, self.arcs[c.arc].name, self.lines[c.line].name, c.distance.abs())); }
2041        for c in &self.distance_aa_ce_ce { out.push(format!("C{}: distance {}.center {}.center = {}", c.nid, self.arcs[c.a].name, self.arcs[c.b].name, c.distance.abs())); }
2042        for c in &self.distance_aa_ce_s { out.push(format!("C{}: distance {}.center {}.start = {}", c.nid, self.arcs[c.a].name, self.arcs[c.b].name, c.distance.abs())); }
2043        for c in &self.distance_aa_ce_e { out.push(format!("C{}: distance {}.center {}.end = {}", c.nid, self.arcs[c.a].name, self.arcs[c.b].name, c.distance.abs())); }
2044        for c in &self.distance_aa_s_ce { out.push(format!("C{}: distance {}.start {}.center = {}", c.nid, self.arcs[c.a].name, self.arcs[c.b].name, c.distance.abs())); }
2045        for c in &self.distance_aa_s_s { out.push(format!("C{}: distance {}.start {}.start = {}", c.nid, self.arcs[c.a].name, self.arcs[c.b].name, c.distance.abs())); }
2046        for c in &self.distance_aa_s_e { out.push(format!("C{}: distance {}.start {}.end = {}", c.nid, self.arcs[c.a].name, self.arcs[c.b].name, c.distance.abs())); }
2047        for c in &self.distance_aa_e_ce { out.push(format!("C{}: distance {}.end {}.center = {}", c.nid, self.arcs[c.a].name, self.arcs[c.b].name, c.distance.abs())); }
2048        for c in &self.distance_aa_e_s { out.push(format!("C{}: distance {}.end {}.start = {}", c.nid, self.arcs[c.a].name, self.arcs[c.b].name, c.distance.abs())); }
2049        for c in &self.distance_aa_e_e { out.push(format!("C{}: distance {}.end {}.end = {}", c.nid, self.arcs[c.a].name, self.arcs[c.b].name, c.distance.abs())); }
2050        for c in &self.distance_concentric {
2051            out.push(format!("C{}: distance {} {} = {} (concentric)", c.nid, self.arcs[c.a].name, self.arcs[c.b].name, c.distance.abs()));
2052        }
2053        // Midpoint variants
2054        for c in &self.midpoint_lp1 { out.push(format!("C{}: midpoint {}.p1 {}", c.nid, self.lines[c.target].name, self.lines[c.line].name)); }
2055        for c in &self.midpoint_lp2 { out.push(format!("C{}: midpoint {}.p2 {}", c.nid, self.lines[c.target].name, self.lines[c.line].name)); }
2056        for c in &self.midpoint_arc_start { out.push(format!("C{}: midpoint {}.start {}", c.nid, self.arcs[c.arc].name, self.lines[c.line].name)); }
2057        for c in &self.midpoint_arc_end { out.push(format!("C{}: midpoint {}.end {}", c.nid, self.arcs[c.arc].name, self.lines[c.line].name)); }
2058        // Axis distance constraints
2059        let axis_label = |h: bool| if h { "hdistance" } else { "vdistance" };
2060        for c in &self.axis_distance_ll11 { out.push(format!("C{}: {} {}.p1 {}.p1 = {}", c.nid, axis_label(c.horizontal), self.lines[c.a].name, self.lines[c.b].name, c.distance.abs())); }
2061        for c in &self.axis_distance_ll12 { out.push(format!("C{}: {} {}.p1 {}.p2 = {}", c.nid, axis_label(c.horizontal), self.lines[c.a].name, self.lines[c.b].name, c.distance.abs())); }
2062        for c in &self.axis_distance_ll21 { out.push(format!("C{}: {} {}.p2 {}.p1 = {}", c.nid, axis_label(c.horizontal), self.lines[c.a].name, self.lines[c.b].name, c.distance.abs())); }
2063        for c in &self.axis_distance_ll22 { out.push(format!("C{}: {} {}.p2 {}.p2 = {}", c.nid, axis_label(c.horizontal), self.lines[c.a].name, self.lines[c.b].name, c.distance.abs())); }
2064        for c in &self.axis_distance_lp1 { out.push(format!("C{}: {} {}.p1 {} = {}", c.nid, axis_label(c.horizontal), self.lines[c.line].name, self.point_display_name(c.point), c.distance.abs())); }
2065        for c in &self.axis_distance_lp2 { out.push(format!("C{}: {} {}.p2 {} = {}", c.nid, axis_label(c.horizontal), self.lines[c.line].name, self.point_display_name(c.point), c.distance.abs())); }
2066        for c in &self.axis_distance_arc_center_p { out.push(format!("C{}: {} {}.center {} = {}", c.nid, axis_label(c.horizontal), self.arcs[c.arc].name, self.point_display_name(c.point), c.distance.abs())); }
2067        for c in &self.axis_distance_arc_start_p { out.push(format!("C{}: {} {}.start {} = {}", c.nid, axis_label(c.horizontal), self.arcs[c.arc].name, self.point_display_name(c.point), c.distance.abs())); }
2068        for c in &self.axis_distance_arc_end_p { out.push(format!("C{}: {} {}.end {} = {}", c.nid, axis_label(c.horizontal), self.arcs[c.arc].name, self.point_display_name(c.point), c.distance.abs())); }
2069        for c in &self.axis_distance_arc_center_l1 { out.push(format!("C{}: {} {}.center {}.p1 = {}", c.nid, axis_label(c.horizontal), self.arcs[c.arc].name, self.lines[c.line].name, c.distance.abs())); }
2070        for c in &self.axis_distance_arc_center_l2 { out.push(format!("C{}: {} {}.center {}.p2 = {}", c.nid, axis_label(c.horizontal), self.arcs[c.arc].name, self.lines[c.line].name, c.distance.abs())); }
2071        for c in &self.axis_distance_arc_start_l1 { out.push(format!("C{}: {} {}.start {}.p1 = {}", c.nid, axis_label(c.horizontal), self.arcs[c.arc].name, self.lines[c.line].name, c.distance.abs())); }
2072        for c in &self.axis_distance_arc_start_l2 { out.push(format!("C{}: {} {}.start {}.p2 = {}", c.nid, axis_label(c.horizontal), self.arcs[c.arc].name, self.lines[c.line].name, c.distance.abs())); }
2073        for c in &self.axis_distance_arc_end_l1 { out.push(format!("C{}: {} {}.end {}.p1 = {}", c.nid, axis_label(c.horizontal), self.arcs[c.arc].name, self.lines[c.line].name, c.distance.abs())); }
2074        for c in &self.axis_distance_arc_end_l2 { out.push(format!("C{}: {} {}.end {}.p2 = {}", c.nid, axis_label(c.horizontal), self.arcs[c.arc].name, self.lines[c.line].name, c.distance.abs())); }
2075        for c in &self.axis_distance_aa_ce_ce { out.push(format!("C{}: {} {}.center {}.center = {}", c.nid, axis_label(c.horizontal), self.arcs[c.a].name, self.arcs[c.b].name, c.distance.abs())); }
2076        for c in &self.axis_distance_aa_ce_s { out.push(format!("C{}: {} {}.center {}.start = {}", c.nid, axis_label(c.horizontal), self.arcs[c.a].name, self.arcs[c.b].name, c.distance.abs())); }
2077        for c in &self.axis_distance_aa_ce_e { out.push(format!("C{}: {} {}.center {}.end = {}", c.nid, axis_label(c.horizontal), self.arcs[c.a].name, self.arcs[c.b].name, c.distance.abs())); }
2078        for c in &self.axis_distance_aa_s_ce { out.push(format!("C{}: {} {}.start {}.center = {}", c.nid, axis_label(c.horizontal), self.arcs[c.a].name, self.arcs[c.b].name, c.distance.abs())); }
2079        for c in &self.axis_distance_aa_s_s { out.push(format!("C{}: {} {}.start {}.start = {}", c.nid, axis_label(c.horizontal), self.arcs[c.a].name, self.arcs[c.b].name, c.distance.abs())); }
2080        for c in &self.axis_distance_aa_s_e { out.push(format!("C{}: {} {}.start {}.end = {}", c.nid, axis_label(c.horizontal), self.arcs[c.a].name, self.arcs[c.b].name, c.distance.abs())); }
2081        for c in &self.axis_distance_aa_e_ce { out.push(format!("C{}: {} {}.end {}.center = {}", c.nid, axis_label(c.horizontal), self.arcs[c.a].name, self.arcs[c.b].name, c.distance.abs())); }
2082        for c in &self.axis_distance_aa_e_s { out.push(format!("C{}: {} {}.end {}.start = {}", c.nid, axis_label(c.horizontal), self.arcs[c.a].name, self.arcs[c.b].name, c.distance.abs())); }
2083        for c in &self.axis_distance_aa_e_e { out.push(format!("C{}: {} {}.end {}.end = {}", c.nid, axis_label(c.horizontal), self.arcs[c.a].name, self.arcs[c.b].name, c.distance.abs())); }
2084        out
2085    }
2086
2087    /// and the constraint is: `measured_property - parsed_expr = 0`.
2088    /// Returns Err if the expression fails to parse or references unknown symbols.
2089    pub fn add_expr_dimension(&mut self, kind: DimensionKind, expr_str: &str,
2090                              offset: vect2d, text_along: f64) -> Result<(), String> {
2091        self.validate_expr(expr_str)?;
2092        let parsed = arael_sym::parse(expr_str).unwrap(); // safe: validate_expr checked parse
2093
2094        let name = format!("d{}", self.next_dimension_id);
2095        self.next_dimension_id += 1;
2096
2097        // Build the measured property expression
2098        let dim = Dimension {
2099            kind, value: 0.0, offset, text_along,
2100            name: name.clone(), expr_str: Some(expr_str.to_string()),
2101            broken: false,
2102            derived: false,
2103            range: None,
2104        };
2105        let measured = dim.measured_symbol(self);
2106        self.dimensions.push(dim);
2107
2108        // Residual: measured - expr = 0
2109        let residual = measured - parsed;
2110        self.add_expr_constraint(residual, format!("{} = {}", name, expr_str));
2111        Ok(())
2112    }
2113
2114    /// Rebuild expression constraints and resolve them so they contribute
2115    /// to Jacobian / cost assembly. Needed before calling calc_jacobian or
2116    /// calc_grad_hessian_dense directly (outside of solve()), e.g. for DOF.
2117    pub fn prepare_expr_constraints(&mut self) {
2118        self.rebuild_expr_constraints();
2119        if !self.expr_constraints.is_empty() {
2120            let mut tmp = Vec::new();
2121            self.serialize64(&mut tmp);
2122            let bag = SymbolBag::build(self);
2123            for ec in &mut self.expr_constraints {
2124                ec.resolve(&bag);
2125            }
2126            self.symbol_bag = Some(bag);
2127        }
2128    }
2129
2130
2131    /// Format an informative error message when Hessian decomposition fails.
2132    fn hessian_error_msg(n: usize, hessian: &[f64], err: &str) -> String {
2133        let nan_count = hessian.iter().filter(|v| v.is_nan()).count();
2134        let inf_count = hessian.iter().filter(|v| v.is_infinite()).count();
2135        format!("DOF computation failed: {} ({}x{} Hessian, {} NaN, {} Inf). \
2136                 This likely indicates a solver bug -- please report it.",
2137                err, n, n, nan_count, inf_count)
2138    }
2139
2140    /// Legacy DOF computation via eigendecomposition of the Hessian J^T J.
2141    ///
2142    /// Kept for comparison and for the `dof eigenvalues` diagnostic. Rank
2143    /// detection from Hessian eigenvalues is numerically unstable at high
2144    /// constraint scales because forming J^T J squares the condition
2145    /// number. For routine DOF counting, prefer [`Sketch::compute_dof`]
2146    /// which uses SVD of the Jacobian directly.
2147    ///
2148    /// Uses nalgebra for n<32, faer for n>=32. Benchmark at n=896 (polygon128):
2149    ///   faer eigenvalues-only:    45ms
2150    ///   faer full eigen:          95ms
2151    ///   nalgebra eigenvalues-only: 110ms
2152    ///   nalgebra full eigen:      220ms
2153    pub fn compute_dof_eigenvalues(&mut self, analyze: bool) -> Result<DofResult, String> {
2154        self.compute_dof_eigenvalues_opt(analyze, true)
2155    }
2156
2157    /// Eigenvalue-based DOF analysis with explicit preconditioning flag.
2158    ///
2159    /// When `preconditioned` is true (the default via
2160    /// [`Self::compute_dof_eigenvalues`]), the Hessian is scaled as
2161    /// `H_N = D^{-1} H D^{-1}` where `D = diag(sqrt(diag(H)))`. This is
2162    /// symmetric Jacobi preconditioning -- it preserves the null-space
2163    /// (and hence rank) exactly, but tames the condition number by
2164    /// folding per-parameter scale differences into the diagonal.
2165    /// `J^T J` would otherwise square an already-wide Jacobian
2166    /// condition number, making eigenvalue-based rank detection fail
2167    /// at large sketch scales. Right eigenvectors are back-transformed
2168    /// to raw parameter space as `v_raw = D^{-1} v_N` (then renormed
2169    /// to unit length) so the reported directions remain physical.
2170    ///
2171    /// When `preconditioned` is false, the eigenvalues and eigenvectors
2172    /// are of the raw Hessian -- useful for debugging residual scaling
2173    /// choices.
2174    pub fn compute_dof_eigenvalues_opt(&mut self, analyze: bool, preconditioned: bool) -> Result<DofResult, String> {
2175        // Strip range barriers before the Hessian is assembled -- their
2176        // contribution is state-dependent (zero inside the feasible
2177        // band, non-zero outside) and would make the reported DOF
2178        // swing as geometry crosses a bound. See the matching comment
2179        // in compute_dof. Restored after computation regardless of
2180        // outcome.
2181        self.prepare_expr_constraints();
2182        let saved_ranges: Vec<_> = {
2183            let (kept, ranges): (Vec<_>, Vec<_>) = std::mem::take(&mut self.expr_constraints)
2184                .into_iter().partition(|ec| ec.label != "range");
2185            self.expr_constraints = kept;
2186            ranges
2187        };
2188        let result = self.compute_dof_eigenvalues_opt_inner(analyze, preconditioned);
2189        self.expr_constraints.extend(saved_ranges);
2190        result
2191    }
2192
2193    fn compute_dof_eigenvalues_opt_inner(&mut self, analyze: bool, preconditioned: bool) -> Result<DofResult, String> {
2194        use arael::simple_lm::LmProblem;
2195        let t_total = if TIMING_DEBUG { Some(std::time::Instant::now()) } else { None };
2196
2197        self.update_tangent_flags();
2198        self.update_perpendicular_flags();
2199        self.update_line_dir_flags();
2200
2201        let saved_drift = self.drift_isigma;
2202        self.drift_isigma = 0.0;
2203
2204        let mut params = Vec::new();
2205        self.serialize64(&mut params);
2206        let n = params.len();
2207        if n == 0 {
2208            self.drift_isigma = saved_drift;
2209            return Ok(DofResult { dof: 0, param_names: Vec::new(), eigenvalues: Vec::new(), eigenvectors: Vec::new() });
2210        }
2211
2212        let param_names = if analyze {
2213            let bag = SymbolBag::build(self);
2214            let mut names = vec![String::new(); n];
2215            for (name, &idx) in &bag.param_indices {
2216                let i = idx as usize;
2217                if i < n && names[i].is_empty() { names[i] = name.clone(); }
2218            }
2219            names
2220        } else {
2221            Vec::new()
2222        };
2223
2224        let t_hessian = if TIMING_DEBUG { Some(std::time::Instant::now()) } else { None };
2225        let mut grad = vec![0.0f64; n];
2226        let mut hessian = vec![0.0f64; n * n];
2227        self.calc_grad_hessian_dense(&params, &mut grad, &mut hessian);
2228        self.drift_isigma = saved_drift;
2229        let t_hessian = t_hessian.map(|t| t.elapsed());
2230
2231        // Symmetric Jacobi preconditioning: scale by `sqrt(diag(H))`
2232        // which equals the Jacobian's column L2 norms. Preserves
2233        // null-space exactly (zero eigenvalues stay zero) but narrows
2234        // the non-zero spectrum, so rank detection becomes robust to
2235        // per-parameter scale mismatches. Skipped in `raw` mode.
2236        let d: Vec<f64> = if preconditioned {
2237            (0..n).map(|i| hessian[i * n + i].max(0.0).sqrt().max(1e-15)).collect()
2238        } else {
2239            vec![1.0; n]
2240        };
2241        if preconditioned {
2242            for i in 0..n {
2243                for k in 0..n {
2244                    hessian[i * n + k] /= d[i] * d[k];
2245                }
2246            }
2247        }
2248
2249        // Determine rank via spectral gap in the lower portion of the spectrum.
2250        let t_eigen = if TIMING_DEBUG { Some(std::time::Instant::now()) } else { None };
2251        let rank_from_evs = |evs: &[f64]| -> usize {
2252            let mut sorted: Vec<f64> = evs.iter().map(|v| v.abs()).collect();
2253            sorted.sort_by(|a, b| a.partial_cmp(b).unwrap());
2254            let max_ev = sorted.last().copied().unwrap_or(0.0);
2255            let upper_bound = max_ev * 0.01;
2256            // Same zero-floor trick as rank_from_svs: let near-zero
2257            // eigenvalues participate in the gap search instead of
2258            // being silently skipped when `lo < 1e-20`.
2259            let floor = max_ev * 1e-20;
2260            let mut best_gap = 0.0f64;
2261            let mut best_cut = 0;
2262            for i in 0..sorted.len().saturating_sub(1) {
2263                let lo = sorted[i].max(floor);
2264                let hi = sorted[i + 1].max(floor);
2265                if lo > upper_bound { break; }
2266                let gap = hi / lo;
2267                if gap > best_gap {
2268                    best_gap = gap;
2269                    best_cut = i + 1;
2270                }
2271            }
2272            if best_gap < 1e3 {
2273                best_cut = sorted.iter().filter(|&&v| v < 1e-15).count();
2274            }
2275
2276            evs.len() - best_cut
2277        };
2278        // Eigenvectors come out in the normalised parameter space.
2279        // Back-transform via `v_raw = D^{-1} v_N` (then renormalise to
2280        // unit length) so the reported direction is a physical one
2281        // the user can act on. When preconditioning is off, D is all
2282        // ones and this is a no-op.
2283        let unscale_evecs = |mut evs: Vec<Vec<f64>>| -> Vec<Vec<f64>> {
2284            if !preconditioned { return evs; }
2285            for v in evs.iter_mut() {
2286                for (i, vi) in v.iter_mut().enumerate() {
2287                    *vi /= d[i];
2288                }
2289                let norm: f64 = v.iter().map(|x| x * x).sum::<f64>().sqrt();
2290                if norm > 1e-20 {
2291                    for vi in v.iter_mut() { *vi /= norm; }
2292                }
2293            }
2294            evs
2295        };
2296
2297        let (method, result) = if n < 32 && analyze {
2298            let h = nalgebra::DMatrix::from_row_slice(n, n, &hessian);
2299            let eigen = nalgebra::SymmetricEigen::new(h);
2300            let eigenvalues: Vec<f64> = eigen.eigenvalues.iter().cloned().collect();
2301            let rank = rank_from_evs(&eigenvalues);
2302            let dof = n.saturating_sub(rank);
2303            let eigenvectors: Vec<Vec<f64>> = unscale_evecs((0..n)
2304                .map(|col| eigen.eigenvectors.column(col).iter().cloned().collect())
2305                .collect());
2306            ("nalgebra eigen", DofResult { dof, param_names, eigenvalues, eigenvectors })
2307        } else if n < 32 {
2308            let h = nalgebra::DMatrix::from_row_slice(n, n, &hessian);
2309            let evs = h.symmetric_eigenvalues();
2310            let evs_vec: Vec<f64> = evs.iter().cloned().collect();
2311            let rank = rank_from_evs(&evs_vec);
2312            let dof = n.saturating_sub(rank);
2313            ("nalgebra eigenvalues-only", DofResult { dof, param_names: Vec::new(), eigenvalues: Vec::new(), eigenvectors: Vec::new() })
2314        } else if analyze {
2315            let faer_h = faer::Mat::from_fn(n, n, |i, k| hessian[i * n + k]);
2316            match faer_h.self_adjoint_eigen(faer::Side::Lower) {
2317                Ok(eigen) => {
2318                    let s = eigen.S().column_vector();
2319                    let u = eigen.U();
2320                    let eigenvalues: Vec<f64> = (0..n).map(|i| s[i]).collect();
2321                    let rank = rank_from_evs(&eigenvalues);
2322                    let dof = n.saturating_sub(rank);
2323                    let eigenvectors: Vec<Vec<f64>> = unscale_evecs((0..n)
2324                        .map(|col| (0..n).map(|row| u[(row, col)]).collect())
2325                        .collect());
2326                    ("faer eigen", DofResult { dof, param_names, eigenvalues, eigenvectors })
2327                }
2328                Err(e) => return Err(Self::hessian_error_msg(n, &hessian, &format!("{:?}", e))),
2329            }
2330        } else {
2331            let faer_h = faer::Mat::from_fn(n, n, |i, k| hessian[i * n + k]);
2332            match faer_h.self_adjoint_eigenvalues(faer::Side::Lower) {
2333                Ok(evs) => {
2334                    let evs_vec: Vec<f64> = evs.to_vec();
2335                    let rank = rank_from_evs(&evs_vec);
2336                    let dof = n.saturating_sub(rank);
2337                    ("faer eigenvalues-only", DofResult { dof, param_names: Vec::new(), eigenvalues: Vec::new(), eigenvectors: Vec::new() })
2338                }
2339                Err(e) => return Err(Self::hessian_error_msg(n, &hessian, &format!("{:?}", e))),
2340            }
2341        };
2342        if TIMING_DEBUG {
2343            let t_h = t_hessian.unwrap().as_secs_f64() * 1000.0;
2344            let t_e = t_eigen.unwrap().elapsed().as_secs_f64() * 1000.0;
2345            let t_t = t_total.unwrap().elapsed().as_secs_f64() * 1000.0;
2346            eprintln!("[DOF-EIG] n={} analyze={} method={} hessian={:.3}ms eigen={:.3}ms total={:.3}ms dof={}",
2347                n, analyze, method, t_h, t_e, t_t, result.dof);
2348        }
2349        Ok(result)
2350    }
2351
2352    /// Compute degrees of freedom. When `analyze` is true, also returns
2353    /// parameter names, eigenvalues, and eigenvectors for free direction
2354    /// classification. When false, only the DOF count is computed (fast path).
2355    ///
2356    /// Rank is determined from the singular values of the Jacobian J.
2357    /// The returned `eigenvalues` field holds sigma_i^2 (mathematically
2358    /// equivalent to eigenvalues of J^T J) and `eigenvectors` holds the
2359    /// right singular vectors of J (columns of V, equivalent to
2360    /// eigenvectors of J^T J). SVD is used because forming J^T J squares
2361    /// the condition number and destroys rank-detection precision at
2362    /// high constraint scales (observed: scale=10000 misreports DOF).
2363    ///
2364    /// Uses nalgebra SVD for n<32, faer SVD for n>=32.
2365    pub fn compute_dof(&mut self, analyze: bool) -> Result<DofResult, String> {
2366        let t_total = if TIMING_DEBUG { Some(std::time::Instant::now()) } else { None };
2367
2368        self.prepare_expr_constraints();
2369        self.update_tangent_flags();
2370        self.update_perpendicular_flags();
2371        self.update_line_dir_flags();
2372
2373        let saved_drift = self.drift_isigma;
2374        self.drift_isigma = 0.0;
2375
2376        let mut params = Vec::new();
2377        self.serialize64(&mut params);
2378        let n = params.len();
2379        if n == 0 {
2380            self.drift_isigma = saved_drift;
2381            return Ok(DofResult { dof: 0, param_names: Vec::new(), eigenvalues: Vec::new(), eigenvectors: Vec::new() });
2382        }
2383
2384        let param_names = if analyze {
2385            let bag = SymbolBag::build(self);
2386            let mut names = vec![String::new(); n];
2387            for (name, &idx) in &bag.param_indices {
2388                let i = idx as usize;
2389                if i < n && names[i].is_empty() { names[i] = name.clone(); }
2390            }
2391            names
2392        } else {
2393            Vec::new()
2394        };
2395
2396        let t_jac = if TIMING_DEBUG { Some(std::time::Instant::now()) } else { None };
2397        let mut jacobian = self.calc_jacobian(&params);
2398        self.drift_isigma = saved_drift;
2399        // Strip one-sided barrier rows (range dimensions, tagged
2400        // label = "range") before rank detection. Their Jacobian row
2401        // is zero inside the feasible band and non-zero outside, so
2402        // including them would make the reported DOF swing by one as
2403        // the user drags the geometry across a bound -- e.g. a rect
2404        // with interval width 1..999 would show DOF 2 at width 1 and
2405        // DOF 3 at width 2. The geometric DOF of the sketch is the
2406        // same either way; bounds are inequality constraints that
2407        // shouldn't count.
2408        jacobian.rows.retain(|r| r.label != "range");
2409        let m = jacobian.num_residuals();
2410        let t_jac = t_jac.map(|t| t.elapsed());
2411
2412        if m == 0 {
2413            // No residuals: every parameter is free.
2414            let result = DofResult {
2415                dof: n,
2416                param_names,
2417                eigenvalues: vec![0.0; n],
2418                eigenvectors: (0..n).map(|i| {
2419                    let mut v = vec![0.0; n];
2420                    v[i] = 1.0;
2421                    v
2422                }).collect(),
2423            };
2424            self.cached_dof = Some(result.dof);
2425            return Ok(result);
2426        }
2427
2428
2429        let t_svd = if TIMING_DEBUG { Some(std::time::Instant::now()) } else { None };
2430        // Determine rank from the singular value spectrum of the
2431        // current-params Jacobian. The gap algorithm mirrors the old
2432        // eigenvalue version but operates on sigma directly: a sharp
2433        // jump between "near zero" and "constrained" shows up as a
2434        // large gap ratio, without the condition-number squaring
2435        // that killed precision in the J^T J approach.
2436        //
2437        // This is the *instantaneous* rank at the current params,
2438        // which can over-count DOF at configurations where two
2439        // Jacobian rows happen to be tangent-aligned -- classic
2440        // example: `hdistance L0.p2 L1.p2 = X` at a horizontal L1
2441        // has a row parallel to the length row at that instant, so
2442        // adding it doesn't lower the SVD rank even though it's
2443        // structurally a new constraint (rotate the triangle and it
2444        // decouples).
2445        //
2446        // The natural fix is to SVD a slightly-perturbed copy of the
2447        // params instead ("structural rank"); we tried it and backed
2448        // it out. Even tiny perturbations introduce a new spectrum
2449        // of small-but-nonzero sigmas in real sketches (robot.cmd
2450        // produced a sigma at 4e-5 with rel 5e-10 of max), and the
2451        // rank threshold that would classify the horizontal-line
2452        // sigma as "real" also classifies the robot sigma as "real",
2453        // reporting DOF=2 when the sketch actually has DOF=3. We
2454        // couldn't find a single threshold that handled both cases
2455        // without breaking another. A correct fix would track a
2456        // *structural* rank delta per constraint at creation time
2457        // and sum those instead of trusting the instantaneous SVD;
2458        // that's left for future work. The blocker analysis uses
2459        // this same current-params Jacobian so its row-span check
2460        // agrees with the DOF check that triggered it.
2461        let rank_from_svs = |svs: &[f64]| -> usize {
2462            let mut sorted: Vec<f64> = svs.iter().map(|v| v.abs()).collect();
2463            sorted.sort_by(|a, b| a.partial_cmp(b).unwrap());
2464            let max_sv = sorted.last().copied().unwrap_or(0.0);
2465            let upper_bound = max_sv * 0.01;
2466            // Floor near-zero sigmas at `max_sv * 1e-20` so they
2467            // participate in the gap search instead of being skipped.
2468            // No plausible real sigma lives below `max_sv * 1e-15`, so a
2469            // floor at `1e-20` leaves 5 orders of magnitude of headroom
2470            // above the floor for real rank sigmas to remain resolvable
2471            // while letting the zero-to-first-real gap compete on its
2472            // own merits. Multiple adjacent zeros collapse to ratio ~1
2473            // against each other (both pulled to the same floor),
2474            // leaving only the last-zero -> first-real transition as a
2475            // huge-gap candidate -- which is exactly what identifies the
2476            // rank cut.
2477            let floor = max_sv * 1e-20;
2478            let mut best_gap = 0.0f64;
2479            let mut best_cut = 0;
2480            for i in 0..sorted.len().saturating_sub(1) {
2481                let lo = sorted[i].max(floor);
2482                let hi = sorted[i + 1].max(floor);
2483                if lo > upper_bound { break; }
2484                let gap = hi / lo;
2485                if gap > best_gap {
2486                    best_gap = gap;
2487                    best_cut = i + 1;
2488                }
2489            }
2490            if best_gap < 1e3 {
2491                best_cut = sorted.iter().filter(|&&v| v < 1e-15).count();
2492            }
2493            svs.len() - best_cut
2494        };
2495        // Column-normalise the Jacobian before running rank detection.
2496        // Per-parameter scale differences (angular columns with
2497        // `radius * isigma` entries vs length columns with `isigma`
2498        // entries) can span >= 10^10, which makes the raw SVD spectrum
2499        // span the same range -- and the gap algorithm then picks a
2500        // spurious intra-rank gap on ill-conditioned sketches. The
2501        // column-normalised SVD preserves the Jacobian's null-space
2502        // (rank) but has a spectrum bounded by the row-space geometry,
2503        // where the gap algorithm is robust. Raw SVD remains available
2504        // via `dof singular raw` for residual-design debugging.
2505        let rank = rank_from_svs(&jacobian.singular_values_column_normalised());
2506        let dof = n.saturating_sub(rank);
2507
2508        let result = if analyze {
2509            let svd = jacobian.svd();
2510            let svs = &svd.singular_values;
2511            let k = svs.len();
2512            let mut eigenvalues = vec![0.0; n];
2513            let mut eigenvectors: Vec<Vec<f64>> = (0..n).map(|_| vec![0.0; n]).collect();
2514            for i in 0..k {
2515                eigenvalues[i] = svs[i] * svs[i];
2516                // svd.v is n x k row-major (right singular vectors as columns).
2517                // Column i contains the direction for singular value svs[i];
2518                // copy it into eigenvectors[i] (direction-major, param-indexed).
2519                for row in 0..n {
2520                    eigenvectors[i][row] = svd.v[row * k + i];
2521                }
2522            }
2523            // Fill any remaining slots with orthogonal basis vectors (not
2524            // mathematically meaningful beyond rank, but keeps dimensions
2525            // consistent for callers that iterate).
2526            for i in k..n {
2527                eigenvectors[i][i] = 1.0;
2528            }
2529            DofResult { dof, param_names, eigenvalues, eigenvectors }
2530        } else {
2531            DofResult { dof, param_names: Vec::new(), eigenvalues: Vec::new(), eigenvectors: Vec::new() }
2532        };
2533        let method = if n < 32 { "Jacobian::svd (nalgebra)" } else { "Jacobian::svd (faer)" };
2534
2535        if TIMING_DEBUG {
2536            let t_j = t_jac.unwrap().as_secs_f64() * 1000.0;
2537            let t_s = t_svd.unwrap().elapsed().as_secs_f64() * 1000.0;
2538            let t_t = t_total.unwrap().elapsed().as_secs_f64() * 1000.0;
2539            eprintln!("[DOF] m={} n={} analyze={} method={} jacobian={:.3}ms svd={:.3}ms total={:.3}ms dof={}",
2540                m, n, analyze, method, t_j, t_s, t_t, result.dof);
2541        }
2542        self.cached_dof = Some(result.dof);
2543        Ok(result)
2544    }
2545
2546    /// Return cached DOF or compute it (count only, no eigenvector analysis).
2547    pub fn dof(&mut self) -> Result<usize, String> {
2548        if let Some(d) = self.cached_dof { return Ok(d); }
2549        let result = self.compute_dof(false)?;
2550        Ok(result.dof)
2551    }
2552
2553    /// Update tangent_la shared-endpoint flags by scanning coincident collections.
2554    pub fn update_tangent_flags(&mut self) {
2555        for t in &mut self.tangent_la {
2556            t.p1_arc_start = self.coincident_lp1_arc_start.iter().any(|c| c.line == t.line && c.arc == t.arc);
2557            t.p1_arc_end = self.coincident_lp1_arc_end.iter().any(|c| c.line == t.line && c.arc == t.arc);
2558            t.p2_arc_start = self.coincident_lp2_arc_start.iter().any(|c| c.line == t.line && c.arc == t.arc);
2559            t.p2_arc_end = self.coincident_lp2_arc_end.iter().any(|c| c.line == t.line && c.arc == t.arc);
2560            // Compute dir_sign for shared-endpoint directed tangent constraint.
2561            // Only set once (when 0.0) to remember the initial tangent direction;
2562            // recomputing each time would follow perturbations instead of preventing flips.
2563            let shared = t.p1_arc_start || t.p1_arc_end || t.p2_arc_start || t.p2_arc_end;
2564            if shared && t.dir_sign.is_nan() {
2565                let l = &self.lines[t.line];
2566                let a = &self.arcs[t.arc];
2567                // Compute arc endpoint and tangent from arc parameters
2568                let angle = if t.p1_arc_start || t.p2_arc_start { a.start_angle.value } else { a.end_angle.value };
2569                let cr = a.rotation.value.cos();
2570                let sr = a.rotation.value.sin();
2571                let ct = angle.cos();
2572                let st = angle.sin();
2573                let ax = a.center.value.x + a.radius.value * ct * cr - a.radius_b.value * st * sr;
2574                let ay = a.center.value.y + a.radius.value * ct * sr + a.radius_b.value * st * cr;
2575                let tx = -a.radius.value * st * cr - a.radius_b.value * ct * sr;
2576                let ty = -a.radius.value * st * sr + a.radius_b.value * ct * cr;
2577                // Direction from arc endpoint to the line's other end
2578                let (dx, dy) = if t.p1_arc_start || t.p1_arc_end {
2579                    (l.p2.value.x - ax, l.p2.value.y - ay)
2580                } else {
2581                    (l.p1.value.x - ax, l.p1.value.y - ay)
2582                };
2583                let dot = dx * tx + dy * ty;
2584                t.dir_sign = if dot >= 0.0 { 1.0 } else { -1.0 };
2585            }
2586        }
2587    }
2588
2589    /// Compute dir_sign for perpendicular constraints on first use.
2590    pub fn update_perpendicular_flags(&mut self) {
2591        for p in &mut self.perpendicular {
2592            if p.dir_sign.is_nan() {
2593                let la = &self.lines[p.a];
2594                let lb = &self.lines[p.b];
2595                let cross = (la.p2.value.x - la.p1.value.x) * (lb.p2.value.y - lb.p1.value.y)
2596                          - (la.p2.value.y - la.p1.value.y) * (lb.p2.value.x - lb.p1.value.x);
2597                p.dir_sign = if cross >= 0.0 { 1.0 } else { -1.0 };
2598            }
2599        }
2600    }
2601
2602    /// Initialize h_dir_sign / v_dir_sign on lines with horizontal/vertical
2603    /// constraints the first time they're seen. Sketches loaded from JSON
2604    /// that predate the dir_sign fields deserialize with NaN; this backfills
2605    /// them from whatever orientation the loaded geometry has.
2606    pub fn update_line_dir_flags(&mut self) {
2607        for l in self.lines.iter_mut() {
2608            if l.constraints.horizontal && l.constraints.h_dir_sign.is_nan() {
2609                let dx = l.p2.value.x - l.p1.value.x;
2610                l.constraints.h_dir_sign = if dx >= 0.0 { 1.0 } else { -1.0 };
2611            }
2612            if l.constraints.vertical && l.constraints.v_dir_sign.is_nan() {
2613                let dy = l.p2.value.y - l.p1.value.y;
2614                l.constraints.v_dir_sign = if dy >= 0.0 { 1.0 } else { -1.0 };
2615            }
2616        }
2617    }
2618
2619    /// Add hidden helper Points coincident with every "free" Line/Arc
2620    /// endpoint -- one that doesn't already appear in any coincident_*
2621    /// constraint linking it to another entity. Returns a state token
2622    /// that must be passed to `remove_drag_auto_anchors` at drag-end.
2623    ///
2624    /// # Why this exists (a hack we couldn't avoid)
2625    ///
2626    /// In long equal-length chains, the `length L = c` Jacobian is
2627    /// rank-deficient when a segment is perpendicular to the chain's
2628    /// drag direction (`sqrt(dx^2 + dy^2)` has zero derivative w.r.t.
2629    /// `dx` when `dx ~ 0`). At that pose nothing in the constraint
2630    /// system pulls a free endpoint along the chain's translation
2631    /// axis: `length` doesn't (rank-deficient), `equal` doesn't (only
2632    /// couples lengths), `coincident` only acts on joined endpoints,
2633    /// and the `drift` regularizer is six orders of magnitude weaker
2634    /// than the constraints. Result: per-frame drag translates the
2635    /// chain forward but leaves the free endpoint behind, so the
2636    /// last segment rotates 90 degrees away from the chain's axis,
2637    /// then snaps back, then flips again -- visibly unstable.
2638    ///
2639    /// Boosting `drift_isigma` is the "obvious" fix and doesn't work
2640    /// at long chains (the boost needed to lock down 20 articulation
2641    /// modes also makes the chain feel rigid). Adding a free Point
2642    /// coincident with the free endpoint -- the user's empirical
2643    /// workaround -- does fix it: the new Point's drift residual
2644    /// gets folded into the soft chain eigenmode through the
2645    /// coincident, so the Newton step now displaces the free
2646    /// endpoint along with the chain instead of leaving it stationary.
2647    /// We don't fully understand why this works at chain20 but a
2648    /// straight drift boost doesn't (Hessian spectra differ only in
2649    /// the 4th decimal at start-of-solve), only that it does. So we
2650    /// automate the workaround.
2651    ///
2652    /// # Mechanics
2653    ///
2654    /// Walks every coincident_* vec to identify joined endpoints,
2655    /// then for each Line p1/p2 (optimizable, not joined) and each
2656    /// Arc start/end (optimizable, not joined) pushes a hidden helper
2657    /// `Point` plus a matching `coincident_lp1`/`coincident_lp2`/
2658    /// `coincident_arc_start`/`coincident_arc_end`. The helper is
2659    /// placed at the endpoint's current position with a 0.001 offset
2660    /// so the parser doesn't reject it as a no-op coincidence.
2661    ///
2662    /// Call AFTER pushing the drag apparatus onto the sketch so the
2663    /// dragged endpoint is already "joined" (to the drag helper) and
2664    /// won't get a redundant auto-anchor.
2665    pub fn add_drag_auto_anchors(&mut self) -> DragAutoAnchorState {
2666        // Arc endpoint world position: c + R * cos(t) * cos(rot) - R_b * sin(t) * sin(rot),
2667        //                              c + R * cos(t) * sin(rot) + R_b * sin(t) * cos(rot)
2668        // where t is the start_angle or end_angle.
2669        fn arc_pt(a: &Arc, angle: f64) -> vect2d {
2670            let ct = angle.cos();
2671            let st = angle.sin();
2672            let cr = a.rotation.value.cos();
2673            let sr = a.rotation.value.sin();
2674            vect2d::new(
2675                a.center.value.x + a.radius.value * ct * cr - a.radius_b.value * st * sr,
2676                a.center.value.y + a.radius.value * ct * sr + a.radius_b.value * st * cr,
2677            )
2678        }
2679
2680        let mut state = DragAutoAnchorState {
2681            helper_points: std::vec::Vec::new(),
2682            coincident_lp1_len: self.coincident_lp1.len(),
2683            coincident_lp2_len: self.coincident_lp2.len(),
2684            coincident_arc_start_len: self.coincident_arc_start.len(),
2685            coincident_arc_end_len: self.coincident_arc_end.len(),
2686        };
2687
2688        // Build sets of joined endpoints by walking every coincident_*
2689        // that touches a Line endpoint or Arc start/end.
2690        let mut joined_lp1: std::collections::HashSet<Ref<Line>> = std::collections::HashSet::new();
2691        let mut joined_lp2: std::collections::HashSet<Ref<Line>> = std::collections::HashSet::new();
2692        let mut joined_arc_start: std::collections::HashSet<Ref<Arc>> = std::collections::HashSet::new();
2693        let mut joined_arc_end: std::collections::HashSet<Ref<Arc>> = std::collections::HashSet::new();
2694
2695        for c in &self.coincident_lp1 { joined_lp1.insert(c.line); }
2696        for c in &self.coincident_lp2 { joined_lp2.insert(c.line); }
2697        for c in &self.coincident_ll11 { joined_lp1.insert(c.a); joined_lp1.insert(c.b); }
2698        for c in &self.coincident_ll12 { joined_lp1.insert(c.a); joined_lp2.insert(c.b); }
2699        for c in &self.coincident_ll21 { joined_lp2.insert(c.a); joined_lp1.insert(c.b); }
2700        for c in &self.coincident_ll22 { joined_lp2.insert(c.a); joined_lp2.insert(c.b); }
2701        for c in &self.coincident_lp1_arc_center { joined_lp1.insert(c.line); }
2702        for c in &self.coincident_lp2_arc_center { joined_lp2.insert(c.line); }
2703        for c in &self.coincident_lp1_arc_start  { joined_lp1.insert(c.line); joined_arc_start.insert(c.arc); }
2704        for c in &self.coincident_lp2_arc_start  { joined_lp2.insert(c.line); joined_arc_start.insert(c.arc); }
2705        for c in &self.coincident_lp1_arc_end    { joined_lp1.insert(c.line); joined_arc_end.insert(c.arc); }
2706        for c in &self.coincident_lp2_arc_end    { joined_lp2.insert(c.line); joined_arc_end.insert(c.arc); }
2707        for c in &self.coincident_arc_start { joined_arc_start.insert(c.arc); }
2708        for c in &self.coincident_arc_end   { joined_arc_end.insert(c.arc); }
2709        for c in &self.coincident_arc_center_start { joined_arc_start.insert(c.b); }
2710        for c in &self.coincident_arc_start_center { joined_arc_start.insert(c.a); }
2711        for c in &self.coincident_arc_center_end { joined_arc_end.insert(c.b); }
2712        for c in &self.coincident_arc_end_center { joined_arc_end.insert(c.a); }
2713        for c in &self.coincident_arc_start_start { joined_arc_start.insert(c.a); joined_arc_start.insert(c.b); }
2714        for c in &self.coincident_arc_start_end   { joined_arc_start.insert(c.a); joined_arc_end.insert(c.b); }
2715        for c in &self.coincident_arc_end_start   { joined_arc_end.insert(c.a); joined_arc_start.insert(c.b); }
2716        for c in &self.coincident_arc_end_end     { joined_arc_end.insert(c.a); joined_arc_end.insert(c.b); }
2717
2718        // Lines
2719        let line_refs: std::vec::Vec<Ref<Line>> = self.lines.refs().collect();
2720        for r in line_refs {
2721            let l = &self.lines[r];
2722            if l.p1.optimize && !joined_lp1.contains(&r) {
2723                let pos = vect2d::new(l.p1.value.x + 0.001, l.p1.value.y);
2724                let p = self.add_helper_point(pos);
2725                self.coincident_lp1.push(CoincidentLP1 {
2726                    line: r, point: p, nid: 0, cid: 0, hb: arael::model::CrossBlock::new(),
2727                });
2728                state.helper_points.push(p);
2729            }
2730            let l = &self.lines[r];
2731            if l.p2.optimize && !joined_lp2.contains(&r) {
2732                let pos = vect2d::new(l.p2.value.x + 0.001, l.p2.value.y);
2733                let p = self.add_helper_point(pos);
2734                self.coincident_lp2.push(CoincidentLP2 {
2735                    line: r, point: p, nid: 0, cid: 0, hb: arael::model::CrossBlock::new(),
2736                });
2737                state.helper_points.push(p);
2738            }
2739        }
2740
2741        // Arcs (only optimizable start/end angles -- skip closed circles
2742        // where the angles are fixed)
2743        let arc_refs: std::vec::Vec<Ref<Arc>> = self.arcs.refs().collect();
2744        for r in arc_refs {
2745            let a = &self.arcs[r];
2746            if a.start_angle.optimize && !joined_arc_start.contains(&r) {
2747                let pt = arc_pt(a, a.start_angle.value);
2748                let pos = vect2d::new(pt.x + 0.001, pt.y);
2749                let p = self.add_helper_point(pos);
2750                self.coincident_arc_start.push(CoincidentArcStart {
2751                    point: p, arc: r, nid: 0, cid: 0, hb: arael::model::CrossBlock::new(),
2752                });
2753                state.helper_points.push(p);
2754            }
2755            let a = &self.arcs[r];
2756            if a.end_angle.optimize && !joined_arc_end.contains(&r) {
2757                let pt = arc_pt(a, a.end_angle.value);
2758                let pos = vect2d::new(pt.x + 0.001, pt.y);
2759                let p = self.add_helper_point(pos);
2760                self.coincident_arc_end.push(CoincidentArcEnd {
2761                    point: p, arc: r, nid: 0, cid: 0, hb: arael::model::CrossBlock::new(),
2762                });
2763                state.helper_points.push(p);
2764            }
2765        }
2766
2767        state
2768    }
2769
2770    /// Roll back the auto-anchors set up by `add_drag_auto_anchors`.
2771    pub fn remove_drag_auto_anchors(&mut self, state: DragAutoAnchorState) {
2772        self.coincident_lp1.truncate(state.coincident_lp1_len);
2773        self.coincident_lp2.truncate(state.coincident_lp2_len);
2774        self.coincident_arc_start.truncate(state.coincident_arc_start_len);
2775        self.coincident_arc_end.truncate(state.coincident_arc_end_len);
2776        for p in state.helper_points {
2777            self.points.remove(p);
2778        }
2779    }
2780
2781    /// Solve the sketch constraints using Levenberg-Marquardt.
2782    /// Uses sparse faer Cholesky for n >= 64 params, dense Cholesky otherwise.
2783    /// When starting cost is high, uses graduated optimization (1% -> 10% ->
2784    /// 100% constraint strength) to avoid ill-conditioning from the large
2785    /// constraint/drift sigma ratio.
2786    pub fn solve(&mut self) -> arael::simple_lm::LmResult<f64> {
2787        use arael::simple_lm::LmProblem;
2788        // Rebuild expression constraints from dimensions with expr_str
2789        // (needed after load/undo since expr_constraints is not serialized)
2790        self.rebuild_expr_constraints();
2791        self.update_tangent_flags();
2792        self.update_perpendicular_flags();
2793        self.update_line_dir_flags();
2794
2795        let mut params64: std::vec::Vec<f64> = std::vec::Vec::new();
2796        self.serialize64(&mut params64);
2797        let n = params64.len();
2798
2799        if n == 0 {
2800            return arael::simple_lm::LmResult {
2801                x: params64, start_cost: 0.0, end_cost: 0.0, iterations: 0,
2802            };
2803        }
2804
2805        // Build symbol bag and resolve expression constraints
2806        if !self.expr_constraints.is_empty() {
2807            let bag = SymbolBag::build(self);
2808            for ec in &mut self.expr_constraints {
2809                ec.resolve(&bag);
2810            }
2811            self.symbol_bag = Some(bag);
2812        }
2813
2814        // Compute starting cost to decide strategy
2815        let start_cost = self.calc_cost(&params64);
2816
2817        // Graduated optimization: when starting cost is high, the Hessian
2818        // condition number (constraint_isigma/drift_isigma)^2 can make LM
2819        // oscillate. Solve with reduced constraint strength first to get
2820        // close to the solution, then ramp up to full strength.
2821        let full_isigma = self.constraint_isigma;
2822        let graduated = start_cost > n as f64 * 1e-3;
2823        let stages: &[f64] = if graduated {
2824            &[0.01, 0.1, 1.0]
2825        } else {
2826            &[1.0]
2827        };
2828
2829        let mut total_iters = 0usize;
2830        let mut result = arael::simple_lm::LmResult {
2831            x: params64.clone(),
2832            start_cost,
2833            end_cost: start_cost,
2834            iterations: 0,
2835        };
2836
2837        for &scale in stages {
2838            self.constraint_isigma = full_isigma * scale;
2839
2840            let mut params = std::vec::Vec::new();
2841            self.serialize64(&mut params);
2842            let cost = self.calc_cost(&params);
2843
2844            let lambda = if cost > 1.0 {
2845                (cost * 1e-6).clamp(1e-4, 1.0)
2846            } else {
2847                1e-6
2848            };
2849
2850            let config = arael::simple_lm::LmConfig::<f64> {
2851                initial_lambda: lambda,
2852                abs_precision: 1e-6,
2853                rel_precision: 1e-4,
2854                cost_threshold: n as f64 * 1e-6,
2855                min_iters: if cost > (n as f64 * 1e-4) { 32 } else { 8 },
2856                verbose: self.verbose,
2857                ..Default::default()
2858            };
2859            let stage_result = if n >= 64 {
2860                arael::simple_lm::solve_sparse_faer(&params, self, &config)
2861            } else {
2862                arael::simple_lm::solve(&params, self, &config)
2863            };
2864            self.deserialize64(&stage_result.x);
2865            total_iters += stage_result.iterations;
2866            result.end_cost = stage_result.end_cost;
2867            result.x = stage_result.x;
2868        }
2869
2870        self.constraint_isigma = full_isigma;
2871        self.normalise_ellipse_rotations();
2872        self.update_expr_dim_values();
2873        result.iterations = total_iters;
2874        result
2875    }
2876
2877    /// Wrap every ellipse's `rotation` param into `(-pi, pi]`. Angles
2878    /// can accumulate integer multiples of `2*pi` across repeated drags
2879    /// or expression-driven updates; without this post-pass the solver
2880    /// can settle on rotations like `3.5*pi` that are numerically fine
2881    /// but render as nonsense in `info` / `list` output. Wrapping by
2882    /// whole `2*pi` turns leaves every geometric quantity (point
2883    /// positions, tangents, curvatures) unchanged; the drift anchor is
2884    /// the committed `.value` itself, so next solve's drift picks up
2885    /// the wrapped angle as its new anchor automatically.
2886    fn normalise_ellipse_rotations(&mut self) {
2887        use arael::utils::rad2rad;
2888        let refs: std::vec::Vec<_> = self.arcs.refs().collect();
2889        for r in refs {
2890            let a = &mut self.arcs[r];
2891            if !a.is_ellipse { continue; }
2892            // Only `rotation` is normalised; `start_angle` / `end_angle`
2893            // parameterise the traversal of the ellipse and may span more
2894            // than 2*pi (e.g. a "large" elliptic arc where the sweep
2895            // itself exceeds pi). Wrapping them would silently change
2896            // which side of the ellipse is drawn.
2897            a.rotation.value = rad2rad(a.rotation.value);
2898        }
2899    }
2900
2901    /// Evaluate expression/derived dimensions and user params, cache their computed values.
2902    pub fn update_expr_dim_values(&mut self) {
2903        let has_work = self.dimensions.iter().any(|d| d.expr_str.is_some() || d.derived || d.range.is_some())
2904            || self.user_params.iter().any(|p| !p.broken);
2905        if !has_work { return; }
2906        let bag = SymbolBag::build(self);
2907        let mut params = Vec::new();
2908        self.serialize64(&mut params);
2909        let vars = bag.eval_vars(&params);
2910        // Update user params first (dims may reference them)
2911        for p in &mut self.user_params {
2912            if p.broken { continue; }
2913            if p.expr_str.trim().parse::<f64>().is_ok() { continue; }
2914            if let Ok(parsed) = arael_sym::parse(&p.expr_str) {
2915                let expanded = expr_constraint::expand_derived(&parsed, &bag);
2916                match expanded.eval(&vars) {
2917                    Ok(val) => p.value = val,
2918                    Err(_) => p.broken = true,
2919                }
2920            }
2921        }
2922        for dim in &mut self.dimensions {
2923            if dim.broken { continue; }
2924            if let Some(ref expr_str) = dim.expr_str
2925                && let Ok(parsed) = arael_sym::parse(expr_str) {
2926                    let expanded = expr_constraint::expand_derived(&parsed, &bag);
2927                    match expanded.eval(&vars) {
2928                        Ok(val) => dim.value = val,
2929                        Err(_) => dim.broken = true,
2930                    }
2931                }
2932        }
2933        // Update derived numeric dims from measured geometry
2934        let derived_vals: Vec<(usize, f64)> = (0..self.dimensions.len())
2935            .filter(|&i| self.dimensions[i].derived && self.dimensions[i].expr_str.is_none() && !self.dimensions[i].broken)
2936            .filter_map(|i| {
2937                let measured = self.dimensions[i].measured_symbol(self);
2938                let expanded = expr_constraint::expand_derived(&measured, &bag);
2939                expanded.eval(&vars).ok().map(|v| (i, v))
2940            })
2941            .collect();
2942        for (i, val) in derived_vals {
2943            self.dimensions[i].value = val;
2944        }
2945        // Range dimensions track the current measured value in `value` for
2946        // display (the bound itself lives in `range`). Same eval shape as
2947        // the derived-dim update above.
2948        let range_vals: Vec<(usize, f64)> = (0..self.dimensions.len())
2949            .filter(|&i| self.dimensions[i].range.is_some())
2950            .filter_map(|i| {
2951                let measured = self.dimensions[i].measured_symbol(self);
2952                let expanded = expr_constraint::expand_derived(&measured, &bag);
2953                expanded.eval(&vars).ok().map(|v| (i, v))
2954            })
2955            .collect();
2956        for (i, val) in range_vals {
2957            self.dimensions[i].value = val;
2958        }
2959    }
2960}
2961
2962#[cfg(test)]
2963mod jacobian_tests {
2964    use super::*;
2965    use arael::simple_lm::LmProblem;
2966    use arael::vect::vect2d;
2967
2968    /// Build a sketch with lines, coincident constraint, and an expression
2969    /// dimension, then validate Jacobian against Hessian and cost.
2970    fn make_test_sketch() -> (Sketch, Vec<f64>) {
2971        let mut sketch = Sketch::new();
2972        let l0 = sketch.add_line(vect2d::new(0.0, 0.0), vect2d::new(3.0, 0.0));
2973        let l1 = sketch.add_line(vect2d::new(3.0, 0.0), vect2d::new(5.0, 2.0));
2974        // Coincident: L0.p2 == L1.p1
2975        sketch.coincident_ll21.push(CoincidentLL21 {
2976            a: l0,
2977            b: l1,
2978            nid: 0, cid: 0,
2979            hb: arael::model::CrossBlock::new(),
2980        });
2981        // Length dimension on L0 (creates an expression constraint)
2982        sketch.lines[l0].constraints.has_length = true;
2983        sketch.lines[l0].constraints.length = 5.0;
2984        sketch.dimensions.push(Dimension {
2985            kind: DimensionKind::LineLength(l0),
2986            value: 5.0, offset: vect2d::new(0.0, 1.0), text_along: 0.0,
2987            name: "d0".into(), expr_str: None, broken: false, derived: false,
2988            range: None,
2989        });
2990
2991        sketch.prepare_expr_constraints();
2992        let mut params = Vec::new();
2993        sketch.serialize64(&mut params);
2994        (sketch, params)
2995    }
2996
2997    #[test]
2998    fn sketch_jacobian_cost_matches() {
2999        let (mut sketch, mut params) = make_test_sketch();
3000        // Perturb so residuals are non-zero
3001        params[0] += 0.1;
3002        params[1] += 0.2;
3003        params[4] -= 0.3;
3004
3005        let j = sketch.calc_jacobian(&params);
3006        let cost_j: f64 = j.rows.iter().map(|r| r.residual * r.residual).sum();
3007        let cost_c = sketch.calc_cost(&params);
3008        assert!(
3009            (cost_j - cost_c).abs() < 1e-10,
3010            "cost mismatch: jacobian={}, calc_cost={}", cost_j, cost_c
3011        );
3012    }
3013
3014    #[test]
3015    fn sketch_jacobian_jtj_matches_hessian() {
3016        let (mut sketch, mut params) = make_test_sketch();
3017        params[0] += 0.1;
3018        params[1] += 0.2;
3019        params[4] -= 0.3;
3020
3021        let j = sketch.calc_jacobian(&params);
3022        let dense = j.to_dense();
3023        let m = j.num_residuals();
3024        let n = j.num_params;
3025
3026        // J^T * J
3027        let mut jtj = vec![0.0f64; n * n];
3028        for i in 0..n {
3029            for k in 0..n {
3030                let mut s = 0.0;
3031                for r in 0..m { s += dense[r * n + i] * dense[r * n + k]; }
3032                jtj[i * n + k] = s;
3033            }
3034        }
3035
3036        // Hessian = 2 * J^T * J
3037        let mut grad = vec![0.0f64; n];
3038        let mut hessian = vec![0.0f64; n * n];
3039        sketch.calc_grad_hessian_dense(&params, &mut grad, &mut hessian);
3040
3041        for i in 0..n {
3042            for k in 0..n {
3043                let expected = 2.0 * jtj[i * n + k];
3044                let actual = hessian[i * n + k];
3045                assert!(
3046                    (expected - actual).abs() < 1e-8,
3047                    "H[{},{}] mismatch: 2*JtJ={}, H={}", i, k, expected, actual
3048                );
3049            }
3050        }
3051    }
3052
3053    #[test]
3054    fn sketch_jacobian_gradient_matches() {
3055        let (mut sketch, mut params) = make_test_sketch();
3056        params[0] += 0.1;
3057        params[1] += 0.2;
3058
3059        let j = sketch.calc_jacobian(&params);
3060        let n = j.num_params;
3061
3062        let mut grad_j = vec![0.0f64; n];
3063        for row in &j.rows {
3064            for &(idx, d) in &row.entries {
3065                grad_j[idx as usize] += 2.0 * row.residual * d;
3066            }
3067        }
3068
3069        let mut grad = vec![0.0f64; n];
3070        let mut hessian = vec![0.0f64; n * n];
3071        sketch.calc_grad_hessian_dense(&params, &mut grad, &mut hessian);
3072
3073        for i in 0..n {
3074            assert!(
3075                (grad_j[i] - grad[i]).abs() < 1e-8,
3076                "grad[{}] mismatch: J={}, GH={}", i, grad_j[i], grad[i]
3077            );
3078        }
3079    }
3080}
arael_sketch_solver/lib.rs

arael_sketch_solver/
lib.rs
