Another project
1use bone_solver::{
2 ConstraintSystem, LineHandle, NewtonConfig, PointHandle, Residual, SolverError,
3 assemble_jacobian, decompose, evaluate_residuals, jacobian_triplets, residual_norm,
4 solve_newton, solve_newton_decomposed, sparsity_pattern,
5};
6use bone_types::{NewtonDamping, Parameter, ParameterIndex, ResidualIndex, SolverResidual};
7
8fn p(idx: u32) -> ParameterIndex {
9 ParameterIndex::new(idx)
10}
11
12fn point(x: u32, y: u32) -> PointHandle {
13 PointHandle { x: p(x), y: p(y) }
14}
15
16fn parameters(values: &[f64]) -> Vec<Parameter> {
17 values.iter().copied().map(Parameter::new).collect()
18}
19
20#[test]
21fn horizontal_residual_row_and_jacobian_match_hand_derivation() {
22 let system = ConstraintSystem::new(
23 parameters(&[0.0, 1.5, 3.0, 4.0]),
24 vec![Residual::Horizontal(LineHandle {
25 a: point(0, 1),
26 b: point(2, 3),
27 })],
28 );
29 assert_eq!(system.row_count(), 1);
30 let values = evaluate_residuals(&system, &[0.0, 1.5, 3.0, 4.0]);
31 assert!((values[0] - (1.5 - 4.0)).abs() < 1e-15);
32 let triplets = jacobian_triplets(&system, &[0.0, 1.5, 3.0, 4.0]);
33 assert_eq!(triplets.len(), 2);
34 let zero = ResidualIndex::new(0);
35 assert!(triplets.contains(&(zero, p(1), 1.0)));
36 assert!(triplets.contains(&(zero, p(3), -1.0)));
37}
38
39#[test]
40fn pin_solves_in_one_newton_step() {
41 let system = ConstraintSystem::new(
42 parameters(&[3.0]),
43 vec![Residual::Pin {
44 param: p(0),
45 target: 7.0,
46 }],
47 );
48 let Ok(out) = solve_newton(&system, NewtonConfig::DEFAULT) else {
49 panic!("pin is linear, must converge");
50 };
51 assert!((out[0].value() - 7.0).abs() < 1e-9);
52}
53
54#[test]
55fn coincident_point_point_is_two_row_residual() {
56 let system = ConstraintSystem::new(
57 parameters(&[0.0, 0.0, 1.0, 2.0]),
58 vec![Residual::CoincidentPointPoint(point(0, 1), point(2, 3))],
59 );
60 assert_eq!(system.row_count(), 2);
61 let r = evaluate_residuals(&system, &[0.0, 0.0, 1.0, 2.0]);
62 assert_eq!(r.len(), 2);
63 assert!((r[0] - (0.0 - 1.0)).abs() < 1e-15);
64 assert!((r[1] - (0.0 - 2.0)).abs() < 1e-15);
65}
66
67#[test]
68fn sparsity_pattern_drops_duplicate_positions() {
69 let system = ConstraintSystem::new(
70 parameters(&[0.0, 0.0]),
71 vec![Residual::Pin {
72 param: p(0),
73 target: 1.0,
74 }],
75 );
76 let pattern = sparsity_pattern(&system);
77 assert_eq!(pattern.rows(), 1);
78 assert_eq!(pattern.cols(), 2);
79 assert_eq!(pattern.entries().len(), 1);
80}
81
82#[test]
83fn default_damping_regularises_rank_deficient_jacobian() {
84 let system = ConstraintSystem::new(
85 parameters(&[0.0, 3.0, 5.0, 7.0]),
86 vec![Residual::Horizontal(LineHandle {
87 a: point(0, 1),
88 b: point(2, 3),
89 })],
90 );
91 let Ok(out) = solve_newton(&system, NewtonConfig::DEFAULT) else {
92 panic!("default damping must regularise rank-deficient J^T J");
93 };
94 assert!((out[1].value() - out[3].value()).abs() < 1e-9);
95}
96
97#[test]
98fn zero_damping_on_rank_deficient_jacobian_does_not_silently_succeed() {
99 let system = ConstraintSystem::new(
100 parameters(&[0.0, 3.0, 5.0, 7.0]),
101 vec![Residual::Horizontal(LineHandle {
102 a: point(0, 1),
103 b: point(2, 3),
104 })],
105 );
106 let cfg = NewtonConfig {
107 damping: NewtonDamping::new(0.0),
108 ..NewtonConfig::DEFAULT
109 };
110 let Err(err) = solve_newton(&system, cfg) else {
111 panic!("damping=0 on rank-deficient system must surface an error");
112 };
113 assert!(
114 matches!(
115 err,
116 SolverError::InvalidSolutionFound { .. } | SolverError::NoSolutionFound { .. }
117 ),
118 "expected InvalidSolutionFound or NoSolutionFound, got {err:?}"
119 );
120}
121
122#[test]
123fn residual_norm_matches_euclidean_length() {
124 let out = evaluate_residuals(
125 &ConstraintSystem::new(
126 parameters(&[3.0]),
127 vec![Residual::Pin {
128 param: p(0),
129 target: 0.0,
130 }],
131 ),
132 &[3.0],
133 );
134 let norm = residual_norm(&out);
135 assert!((norm.value() - 3.0).abs() < 1e-15);
136 let _ = SolverResidual::new(norm.value());
137}
138
139#[test]
140fn decomposed_singular_index_is_translated_to_global_parameter() {
141 let pins: Vec<Residual> = (0..10)
142 .map(|i| Residual::Pin {
143 param: p(i),
144 target: 0.0,
145 })
146 .collect();
147 let residuals: Vec<Residual> = pins
148 .into_iter()
149 .chain(std::iter::once(Residual::Horizontal(LineHandle {
150 a: point(10, 11),
151 b: point(12, 13),
152 })))
153 .collect();
154 let seed = [0.0; 10]
155 .iter()
156 .copied()
157 .chain([0.0, 5.0, 1.0, 7.0])
158 .collect::<Vec<_>>();
159 let system = ConstraintSystem::new(parameters(&seed), residuals);
160 let decomp = decompose(&system);
161 let cfg = NewtonConfig {
162 damping: NewtonDamping::new(0.0),
163 ..NewtonConfig::DEFAULT
164 };
165 match solve_newton_decomposed(&system, &decomp, cfg) {
166 Err(SolverError::InvalidSolutionFound { at }) => {
167 assert!(
168 at.as_usize() >= 10,
169 "expected global parameter index (>=10), got {at:?}"
170 );
171 }
172 Err(SolverError::NoSolutionFound { .. }) => {}
173 other => panic!("expected InvalidSolutionFound or NoSolutionFound, got {other:?}"),
174 }
175}
176
177#[test]
178fn symmetric_residual_is_zero_when_b_is_reflection_of_a() {
179 let system = ConstraintSystem::new(
180 parameters(&[0.0, 1.0, 0.0, -1.0, -2.0, 0.0, 2.0, 0.0]),
181 vec![Residual::Symmetric {
182 a: point(0, 1),
183 b: point(2, 3),
184 axis: LineHandle {
185 a: point(4, 5),
186 b: point(6, 7),
187 },
188 }],
189 );
190 assert_eq!(system.row_count(), 2);
191 let values = evaluate_residuals(&system, &[0.0, 1.0, 0.0, -1.0, -2.0, 0.0, 2.0, 0.0]);
192 assert!(values[0].abs() < 1e-15, "perpendicularity: {}", values[0]);
193 assert!(values[1].abs() < 1e-15, "midpoint-on-axis: {}", values[1]);
194}
195
196#[test]
197fn symmetric_residual_is_nonzero_when_b_drifts_off_reflection() {
198 let system = ConstraintSystem::new(
199 parameters(&[0.0, 1.0, 0.5, -1.0, -2.0, 0.0, 2.0, 0.0]),
200 vec![Residual::Symmetric {
201 a: point(0, 1),
202 b: point(2, 3),
203 axis: LineHandle {
204 a: point(4, 5),
205 b: point(6, 7),
206 },
207 }],
208 );
209 let values = evaluate_residuals(&system, &[0.0, 1.0, 0.5, -1.0, -2.0, 0.0, 2.0, 0.0]);
210 assert!(
211 values[0].abs() > 0.0 || values[1].abs() > 0.0,
212 "drift must surface in at least one residual",
213 );
214}
215
216#[test]
217fn symmetric_solves_to_b_as_reflection_of_a_about_x_axis() {
218 let system = ConstraintSystem::new(
219 parameters(&[0.0, 1.0, 0.5, 1.5, -2.0, 0.0, 2.0, 0.0]),
220 vec![
221 Residual::Pin {
222 param: p(0),
223 target: 0.0,
224 },
225 Residual::Pin {
226 param: p(1),
227 target: 1.0,
228 },
229 Residual::Pin {
230 param: p(4),
231 target: -2.0,
232 },
233 Residual::Pin {
234 param: p(5),
235 target: 0.0,
236 },
237 Residual::Pin {
238 param: p(6),
239 target: 2.0,
240 },
241 Residual::Pin {
242 param: p(7),
243 target: 0.0,
244 },
245 Residual::Symmetric {
246 a: point(0, 1),
247 b: point(2, 3),
248 axis: LineHandle {
249 a: point(4, 5),
250 b: point(6, 7),
251 },
252 },
253 ],
254 );
255 let Ok(out) = solve_newton(&system, NewtonConfig::DEFAULT) else {
256 panic!("symmetric over fixed axis must converge");
257 };
258 assert!(
259 (out[2].value() - 0.0).abs() < 1e-6,
260 "b.x = a.x for reflection over x-axis: got {}",
261 out[2].value()
262 );
263 assert!(
264 (out[3].value() - -1.0).abs() < 1e-6,
265 "b.y = -a.y for reflection over x-axis: got {}",
266 out[3].value()
267 );
268}
269
270#[test]
271fn assemble_jacobian_and_triplets_agree_on_dense_shape() {
272 let system = ConstraintSystem::new(
273 parameters(&[0.0, 0.0, 2.0, 0.0]),
274 vec![Residual::Horizontal(LineHandle {
275 a: point(0, 1),
276 b: point(2, 3),
277 })],
278 );
279 let params = vec![0.0, 0.0, 2.0, 0.0];
280 let dense = assemble_jacobian(&system, ¶ms);
281 assert_eq!(dense.nrows(), 1);
282 assert_eq!(dense.ncols(), 4);
283 let triplets = jacobian_triplets(&system, ¶ms);
284 let sum: f64 = triplets
285 .iter()
286 .map(|(r, c, v)| dense[(r.as_usize(), c.as_usize())] - v)
287 .sum();
288 assert!(sum.abs() < 1e-15);
289}