1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.apache.commons.math.estimation;
19
20 import java.util.ArrayList;
21 import java.util.HashSet;
22
23 import junit.framework.Test;
24 import junit.framework.TestCase;
25 import junit.framework.TestSuite;
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89 @Deprecated
90 public class LevenbergMarquardtEstimatorTest
91 extends TestCase {
92
93 public LevenbergMarquardtEstimatorTest(String name) {
94 super(name);
95 }
96
97 public void testTrivial() throws EstimationException {
98 LinearProblem problem =
99 new LinearProblem(new LinearMeasurement[] {
100 new LinearMeasurement(new double[] {2},
101 new EstimatedParameter[] {
102 new EstimatedParameter("p0", 0)
103 }, 3.0)
104 });
105 LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
106 estimator.estimate(problem);
107 assertEquals(0, estimator.getRMS(problem), 1.0e-10);
108 try {
109 estimator.guessParametersErrors(problem);
110 fail("an exception should have been thrown");
111 } catch (EstimationException ee) {
112
113 } catch (Exception e) {
114 fail("wrong exception caught");
115 }
116 assertEquals(1.5,
117 problem.getUnboundParameters()[0].getEstimate(),
118 1.0e-10);
119 }
120
121 public void testQRColumnsPermutation() throws EstimationException {
122
123 EstimatedParameter[] x = {
124 new EstimatedParameter("p0", 0), new EstimatedParameter("p1", 0)
125 };
126 LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
127 new LinearMeasurement(new double[] { 1.0, -1.0 },
128 new EstimatedParameter[] { x[0], x[1] },
129 4.0),
130 new LinearMeasurement(new double[] { 2.0 },
131 new EstimatedParameter[] { x[1] },
132 6.0),
133 new LinearMeasurement(new double[] { 1.0, -2.0 },
134 new EstimatedParameter[] { x[0], x[1] },
135 1.0)
136 });
137
138 LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
139 estimator.estimate(problem);
140 assertEquals(0, estimator.getRMS(problem), 1.0e-10);
141 assertEquals(7.0, x[0].getEstimate(), 1.0e-10);
142 assertEquals(3.0, x[1].getEstimate(), 1.0e-10);
143
144 }
145
146 public void testNoDependency() throws EstimationException {
147 EstimatedParameter[] p = new EstimatedParameter[] {
148 new EstimatedParameter("p0", 0),
149 new EstimatedParameter("p1", 0),
150 new EstimatedParameter("p2", 0),
151 new EstimatedParameter("p3", 0),
152 new EstimatedParameter("p4", 0),
153 new EstimatedParameter("p5", 0)
154 };
155 LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
156 new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[0] }, 0.0),
157 new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[1] }, 1.1),
158 new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[2] }, 2.2),
159 new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[3] }, 3.3),
160 new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[4] }, 4.4),
161 new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[5] }, 5.5)
162 });
163 LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
164 estimator.estimate(problem);
165 assertEquals(0, estimator.getRMS(problem), 1.0e-10);
166 for (int i = 0; i < p.length; ++i) {
167 assertEquals(0.55 * i, p[i].getEstimate(), 1.0e-10);
168 }
169 }
170
171 public void testOneSet() throws EstimationException {
172
173 EstimatedParameter[] p = {
174 new EstimatedParameter("p0", 0),
175 new EstimatedParameter("p1", 0),
176 new EstimatedParameter("p2", 0)
177 };
178 LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
179 new LinearMeasurement(new double[] { 1.0 },
180 new EstimatedParameter[] { p[0] },
181 1.0),
182 new LinearMeasurement(new double[] { -1.0, 1.0 },
183 new EstimatedParameter[] { p[0], p[1] },
184 1.0),
185 new LinearMeasurement(new double[] { -1.0, 1.0 },
186 new EstimatedParameter[] { p[1], p[2] },
187 1.0)
188 });
189
190 LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
191 estimator.estimate(problem);
192 assertEquals(0, estimator.getRMS(problem), 1.0e-10);
193 assertEquals(1.0, p[0].getEstimate(), 1.0e-10);
194 assertEquals(2.0, p[1].getEstimate(), 1.0e-10);
195 assertEquals(3.0, p[2].getEstimate(), 1.0e-10);
196
197 }
198
199 public void testTwoSets() throws EstimationException {
200 EstimatedParameter[] p = {
201 new EstimatedParameter("p0", 0),
202 new EstimatedParameter("p1", 1),
203 new EstimatedParameter("p2", 2),
204 new EstimatedParameter("p3", 3),
205 new EstimatedParameter("p4", 4),
206 new EstimatedParameter("p5", 5)
207 };
208
209 double epsilon = 1.0e-7;
210 LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
211
212
213 new LinearMeasurement(new double[] { 2.0, 1.0, 4.0 },
214 new EstimatedParameter[] { p[0], p[1], p[3] },
215 2.0),
216 new LinearMeasurement(new double[] { -4.0, -2.0, 3.0, -7.0 },
217 new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
218 -9.0),
219 new LinearMeasurement(new double[] { 4.0, 1.0, -2.0, 8.0 },
220 new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
221 2.0),
222 new LinearMeasurement(new double[] { -3.0, -12.0, -1.0 },
223 new EstimatedParameter[] { p[1], p[2], p[3] },
224 2.0),
225
226
227 new LinearMeasurement(new double[] { epsilon, 1.0 },
228 new EstimatedParameter[] { p[4], p[5] },
229 1.0 + epsilon * epsilon),
230 new LinearMeasurement(new double[] { 1.0, 1.0 },
231 new EstimatedParameter[] { p[4], p[5] },
232 2.0)
233
234 });
235
236 LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
237 estimator.estimate(problem);
238 assertEquals(0, estimator.getRMS(problem), 1.0e-10);
239 assertEquals( 3.0, p[0].getEstimate(), 1.0e-10);
240 assertEquals( 4.0, p[1].getEstimate(), 1.0e-10);
241 assertEquals(-1.0, p[2].getEstimate(), 1.0e-10);
242 assertEquals(-2.0, p[3].getEstimate(), 1.0e-10);
243 assertEquals( 1.0 + epsilon, p[4].getEstimate(), 1.0e-10);
244 assertEquals( 1.0 - epsilon, p[5].getEstimate(), 1.0e-10);
245
246 }
247
248 public void testNonInversible() throws EstimationException {
249
250 EstimatedParameter[] p = {
251 new EstimatedParameter("p0", 0),
252 new EstimatedParameter("p1", 0),
253 new EstimatedParameter("p2", 0)
254 };
255 LinearMeasurement[] m = new LinearMeasurement[] {
256 new LinearMeasurement(new double[] { 1.0, 2.0, -3.0 },
257 new EstimatedParameter[] { p[0], p[1], p[2] },
258 1.0),
259 new LinearMeasurement(new double[] { 2.0, 1.0, 3.0 },
260 new EstimatedParameter[] { p[0], p[1], p[2] },
261 1.0),
262 new LinearMeasurement(new double[] { -3.0, -9.0 },
263 new EstimatedParameter[] { p[0], p[2] },
264 1.0)
265 };
266 LinearProblem problem = new LinearProblem(m);
267
268 LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
269 double initialCost = estimator.getRMS(problem);
270 estimator.estimate(problem);
271 assertTrue(estimator.getRMS(problem) < initialCost);
272 assertTrue(Math.sqrt(m.length) * estimator.getRMS(problem) > 0.6);
273 try {
274 estimator.getCovariances(problem);
275 fail("an exception should have been thrown");
276 } catch (EstimationException ee) {
277
278 } catch (Exception e) {
279 fail("wrong exception caught");
280 }
281 double dJ0 = 2 * (m[0].getResidual() * m[0].getPartial(p[0])
282 + m[1].getResidual() * m[1].getPartial(p[0])
283 + m[2].getResidual() * m[2].getPartial(p[0]));
284 double dJ1 = 2 * (m[0].getResidual() * m[0].getPartial(p[1])
285 + m[1].getResidual() * m[1].getPartial(p[1]));
286 double dJ2 = 2 * (m[0].getResidual() * m[0].getPartial(p[2])
287 + m[1].getResidual() * m[1].getPartial(p[2])
288 + m[2].getResidual() * m[2].getPartial(p[2]));
289 assertEquals(0, dJ0, 1.0e-10);
290 assertEquals(0, dJ1, 1.0e-10);
291 assertEquals(0, dJ2, 1.0e-10);
292
293 }
294
295 public void testIllConditioned() throws EstimationException {
296 EstimatedParameter[] p = {
297 new EstimatedParameter("p0", 0),
298 new EstimatedParameter("p1", 1),
299 new EstimatedParameter("p2", 2),
300 new EstimatedParameter("p3", 3)
301 };
302
303 LinearProblem problem1 = new LinearProblem(new LinearMeasurement[] {
304 new LinearMeasurement(new double[] { 10.0, 7.0, 8.0, 7.0 },
305 new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
306 32.0),
307 new LinearMeasurement(new double[] { 7.0, 5.0, 6.0, 5.0 },
308 new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
309 23.0),
310 new LinearMeasurement(new double[] { 8.0, 6.0, 10.0, 9.0 },
311 new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
312 33.0),
313 new LinearMeasurement(new double[] { 7.0, 5.0, 9.0, 10.0 },
314 new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
315 31.0)
316 });
317 LevenbergMarquardtEstimator estimator1 = new LevenbergMarquardtEstimator();
318 estimator1.estimate(problem1);
319 assertEquals(0, estimator1.getRMS(problem1), 1.0e-10);
320 assertEquals(1.0, p[0].getEstimate(), 1.0e-10);
321 assertEquals(1.0, p[1].getEstimate(), 1.0e-10);
322 assertEquals(1.0, p[2].getEstimate(), 1.0e-10);
323 assertEquals(1.0, p[3].getEstimate(), 1.0e-10);
324
325 LinearProblem problem2 = new LinearProblem(new LinearMeasurement[] {
326 new LinearMeasurement(new double[] { 10.0, 7.0, 8.1, 7.2 },
327 new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
328 32.0),
329 new LinearMeasurement(new double[] { 7.08, 5.04, 6.0, 5.0 },
330 new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
331 23.0),
332 new LinearMeasurement(new double[] { 8.0, 5.98, 9.89, 9.0 },
333 new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
334 33.0),
335 new LinearMeasurement(new double[] { 6.99, 4.99, 9.0, 9.98 },
336 new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
337 31.0)
338 });
339 LevenbergMarquardtEstimator estimator2 = new LevenbergMarquardtEstimator();
340 estimator2.estimate(problem2);
341 assertEquals(0, estimator2.getRMS(problem2), 1.0e-10);
342 assertEquals(-81.0, p[0].getEstimate(), 1.0e-8);
343 assertEquals(137.0, p[1].getEstimate(), 1.0e-8);
344 assertEquals(-34.0, p[2].getEstimate(), 1.0e-8);
345 assertEquals( 22.0, p[3].getEstimate(), 1.0e-8);
346
347 }
348
349 public void testMoreEstimatedParametersSimple() throws EstimationException {
350
351 EstimatedParameter[] p = {
352 new EstimatedParameter("p0", 7),
353 new EstimatedParameter("p1", 6),
354 new EstimatedParameter("p2", 5),
355 new EstimatedParameter("p3", 4)
356 };
357 LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
358 new LinearMeasurement(new double[] { 3.0, 2.0 },
359 new EstimatedParameter[] { p[0], p[1] },
360 7.0),
361 new LinearMeasurement(new double[] { 1.0, -1.0, 1.0 },
362 new EstimatedParameter[] { p[1], p[2], p[3] },
363 3.0),
364 new LinearMeasurement(new double[] { 2.0, 1.0 },
365 new EstimatedParameter[] { p[0], p[2] },
366 5.0)
367 });
368
369 LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
370 estimator.estimate(problem);
371 assertEquals(0, estimator.getRMS(problem), 1.0e-10);
372
373 }
374
375 public void testMoreEstimatedParametersUnsorted() throws EstimationException {
376 EstimatedParameter[] p = {
377 new EstimatedParameter("p0", 2),
378 new EstimatedParameter("p1", 2),
379 new EstimatedParameter("p2", 2),
380 new EstimatedParameter("p3", 2),
381 new EstimatedParameter("p4", 2),
382 new EstimatedParameter("p5", 2)
383 };
384 LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
385 new LinearMeasurement(new double[] { 1.0, 1.0 },
386 new EstimatedParameter[] { p[0], p[1] },
387 3.0),
388 new LinearMeasurement(new double[] { 1.0, 1.0, 1.0 },
389 new EstimatedParameter[] { p[2], p[3], p[4] },
390 12.0),
391 new LinearMeasurement(new double[] { 1.0, -1.0 },
392 new EstimatedParameter[] { p[4], p[5] },
393 -1.0),
394 new LinearMeasurement(new double[] { 1.0, -1.0, 1.0 },
395 new EstimatedParameter[] { p[3], p[2], p[5] },
396 7.0),
397 new LinearMeasurement(new double[] { 1.0, -1.0 },
398 new EstimatedParameter[] { p[4], p[3] },
399 1.0)
400 });
401
402 LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
403 estimator.estimate(problem);
404 assertEquals(0, estimator.getRMS(problem), 1.0e-10);
405 assertEquals(3.0, p[2].getEstimate(), 1.0e-10);
406 assertEquals(4.0, p[3].getEstimate(), 1.0e-10);
407 assertEquals(5.0, p[4].getEstimate(), 1.0e-10);
408 assertEquals(6.0, p[5].getEstimate(), 1.0e-10);
409
410 }
411
412 public void testRedundantEquations() throws EstimationException {
413 EstimatedParameter[] p = {
414 new EstimatedParameter("p0", 1),
415 new EstimatedParameter("p1", 1)
416 };
417 LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
418 new LinearMeasurement(new double[] { 1.0, 1.0 },
419 new EstimatedParameter[] { p[0], p[1] },
420 3.0),
421 new LinearMeasurement(new double[] { 1.0, -1.0 },
422 new EstimatedParameter[] { p[0], p[1] },
423 1.0),
424 new LinearMeasurement(new double[] { 1.0, 3.0 },
425 new EstimatedParameter[] { p[0], p[1] },
426 5.0)
427 });
428
429 LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
430 estimator.estimate(problem);
431 assertEquals(0, estimator.getRMS(problem), 1.0e-10);
432 assertEquals(2.0, p[0].getEstimate(), 1.0e-10);
433 assertEquals(1.0, p[1].getEstimate(), 1.0e-10);
434
435 }
436
437 public void testInconsistentEquations() throws EstimationException {
438 EstimatedParameter[] p = {
439 new EstimatedParameter("p0", 1),
440 new EstimatedParameter("p1", 1)
441 };
442 LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
443 new LinearMeasurement(new double[] { 1.0, 1.0 },
444 new EstimatedParameter[] { p[0], p[1] },
445 3.0),
446 new LinearMeasurement(new double[] { 1.0, -1.0 },
447 new EstimatedParameter[] { p[0], p[1] },
448 1.0),
449 new LinearMeasurement(new double[] { 1.0, 3.0 },
450 new EstimatedParameter[] { p[0], p[1] },
451 4.0)
452 });
453
454 LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
455 estimator.estimate(problem);
456 assertTrue(estimator.getRMS(problem) > 0.1);
457
458 }
459
460 public void testControlParameters() {
461 Circle circle = new Circle(98.680, 47.345);
462 circle.addPoint( 30.0, 68.0);
463 circle.addPoint( 50.0, -6.0);
464 circle.addPoint(110.0, -20.0);
465 circle.addPoint( 35.0, 15.0);
466 circle.addPoint( 45.0, 97.0);
467 checkEstimate(circle, 0.1, 10, 1.0e-14, 1.0e-16, 1.0e-10, false);
468 checkEstimate(circle, 0.1, 10, 1.0e-15, 1.0e-17, 1.0e-10, true);
469 checkEstimate(circle, 0.1, 5, 1.0e-15, 1.0e-16, 1.0e-10, true);
470 circle.addPoint(300, -300);
471 checkEstimate(circle, 0.1, 20, 1.0e-18, 1.0e-16, 1.0e-10, true);
472 }
473
474 private void checkEstimate(EstimationProblem problem,
475 double initialStepBoundFactor, int maxCostEval,
476 double costRelativeTolerance, double parRelativeTolerance,
477 double orthoTolerance, boolean shouldFail) {
478 try {
479 LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
480 estimator.setInitialStepBoundFactor(initialStepBoundFactor);
481 estimator.setMaxCostEval(maxCostEval);
482 estimator.setCostRelativeTolerance(costRelativeTolerance);
483 estimator.setParRelativeTolerance(parRelativeTolerance);
484 estimator.setOrthoTolerance(orthoTolerance);
485 estimator.estimate(problem);
486 assertTrue(! shouldFail);
487 } catch (EstimationException ee) {
488 assertTrue(shouldFail);
489 } catch (Exception e) {
490 fail("wrong exception type caught");
491 }
492 }
493
494 public void testCircleFitting() throws EstimationException {
495 Circle circle = new Circle(98.680, 47.345);
496 circle.addPoint( 30.0, 68.0);
497 circle.addPoint( 50.0, -6.0);
498 circle.addPoint(110.0, -20.0);
499 circle.addPoint( 35.0, 15.0);
500 circle.addPoint( 45.0, 97.0);
501 LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
502 estimator.estimate(circle);
503 assertTrue(estimator.getCostEvaluations() < 10);
504 assertTrue(estimator.getJacobianEvaluations() < 10);
505 double rms = estimator.getRMS(circle);
506 assertEquals(1.768262623567235, Math.sqrt(circle.getM()) * rms, 1.0e-10);
507 assertEquals(69.96016176931406, circle.getRadius(), 1.0e-10);
508 assertEquals(96.07590211815305, circle.getX(), 1.0e-10);
509 assertEquals(48.13516790438953, circle.getY(), 1.0e-10);
510 double[][] cov = estimator.getCovariances(circle);
511 assertEquals(1.839, cov[0][0], 0.001);
512 assertEquals(0.731, cov[0][1], 0.001);
513 assertEquals(cov[0][1], cov[1][0], 1.0e-14);
514 assertEquals(0.786, cov[1][1], 0.001);
515 double[] errors = estimator.guessParametersErrors(circle);
516 assertEquals(1.384, errors[0], 0.001);
517 assertEquals(0.905, errors[1], 0.001);
518
519
520 double cx = circle.getX();
521 double cy = circle.getY();
522 double r = circle.getRadius();
523 for (double d= 0; d < 2 * Math.PI; d += 0.01) {
524 circle.addPoint(cx + r * Math.cos(d), cy + r * Math.sin(d));
525 }
526 estimator = new LevenbergMarquardtEstimator();
527 estimator.estimate(circle);
528 cov = estimator.getCovariances(circle);
529 assertEquals(0.004, cov[0][0], 0.001);
530 assertEquals(6.40e-7, cov[0][1], 1.0e-9);
531 assertEquals(cov[0][1], cov[1][0], 1.0e-14);
532 assertEquals(0.003, cov[1][1], 0.001);
533 errors = estimator.guessParametersErrors(circle);
534 assertEquals(0.004, errors[0], 0.001);
535 assertEquals(0.004, errors[1], 0.001);
536
537 }
538
539 public void testCircleFittingBadInit() throws EstimationException {
540 Circle circle = new Circle(-12, -12);
541 double[][] points = new double[][] {
542 {-0.312967, 0.072366}, {-0.339248, 0.132965}, {-0.379780, 0.202724},
543 {-0.390426, 0.260487}, {-0.361212, 0.328325}, {-0.346039, 0.392619},
544 {-0.280579, 0.444306}, {-0.216035, 0.470009}, {-0.149127, 0.493832},
545 {-0.075133, 0.483271}, {-0.007759, 0.452680}, { 0.060071, 0.410235},
546 { 0.103037, 0.341076}, { 0.118438, 0.273884}, { 0.131293, 0.192201},
547 { 0.115869, 0.129797}, { 0.072223, 0.058396}, { 0.022884, 0.000718},
548 {-0.053355, -0.020405}, {-0.123584, -0.032451}, {-0.216248, -0.032862},
549 {-0.278592, -0.005008}, {-0.337655, 0.056658}, {-0.385899, 0.112526},
550 {-0.405517, 0.186957}, {-0.415374, 0.262071}, {-0.387482, 0.343398},
551 {-0.347322, 0.397943}, {-0.287623, 0.458425}, {-0.223502, 0.475513},
552 {-0.135352, 0.478186}, {-0.061221, 0.483371}, { 0.003711, 0.422737},
553 { 0.065054, 0.375830}, { 0.108108, 0.297099}, { 0.123882, 0.222850},
554 { 0.117729, 0.134382}, { 0.085195, 0.056820}, { 0.029800, -0.019138},
555 {-0.027520, -0.072374}, {-0.102268, -0.091555}, {-0.200299, -0.106578},
556 {-0.292731, -0.091473}, {-0.356288, -0.051108}, {-0.420561, 0.014926},
557 {-0.471036, 0.074716}, {-0.488638, 0.182508}, {-0.485990, 0.254068},
558 {-0.463943, 0.338438}, {-0.406453, 0.404704}, {-0.334287, 0.466119},
559 {-0.254244, 0.503188}, {-0.161548, 0.495769}, {-0.075733, 0.495560},
560 { 0.001375, 0.434937}, { 0.082787, 0.385806}, { 0.115490, 0.323807},
561 { 0.141089, 0.223450}, { 0.138693, 0.131703}, { 0.126415, 0.049174},
562 { 0.066518, -0.010217}, {-0.005184, -0.070647}, {-0.080985, -0.103635},
563 {-0.177377, -0.116887}, {-0.260628, -0.100258}, {-0.335756, -0.056251},
564 {-0.405195, -0.000895}, {-0.444937, 0.085456}, {-0.484357, 0.175597},
565 {-0.472453, 0.248681}, {-0.438580, 0.347463}, {-0.402304, 0.422428},
566 {-0.326777, 0.479438}, {-0.247797, 0.505581}, {-0.152676, 0.519380},
567 {-0.071754, 0.516264}, { 0.015942, 0.472802}, { 0.076608, 0.419077},
568 { 0.127673, 0.330264}, { 0.159951, 0.262150}, { 0.153530, 0.172681},
569 { 0.140653, 0.089229}, { 0.078666, 0.024981}, { 0.023807, -0.037022},
570 {-0.048837, -0.077056}, {-0.127729, -0.075338}, {-0.221271, -0.067526}
571 };
572 for (int i = 0; i < points.length; ++i) {
573 circle.addPoint(points[i][0], points[i][1]);
574 }
575 LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
576 estimator.estimate(circle);
577 assertTrue(estimator.getCostEvaluations() < 15);
578 assertTrue(estimator.getJacobianEvaluations() < 10);
579 assertEquals( 0.030184491196225207, estimator.getRMS(circle), 1.0e-9);
580 assertEquals( 0.2922350065939634, circle.getRadius(), 1.0e-9);
581 assertEquals(-0.15173845023862165, circle.getX(), 1.0e-8);
582 assertEquals( 0.20750021499570379, circle.getY(), 1.0e-8);
583 }
584
585 public void testMath199() {
586 try {
587 QuadraticProblem problem = new QuadraticProblem();
588 problem.addPoint (0, -3.182591015485607, 0.0);
589 problem.addPoint (1, -2.5581184967730577, 4.4E-323);
590 problem.addPoint (2, -2.1488478161387325, 1.0);
591 problem.addPoint (3, -1.9122489313410047, 4.4E-323);
592 problem.addPoint (4, 1.7785661310051026, 0.0);
593 new LevenbergMarquardtEstimator().estimate(problem);
594 fail("an exception should have been thrown");
595 } catch (EstimationException ee) {
596
597 }
598
599 }
600
601 private static class LinearProblem implements EstimationProblem {
602
603 public LinearProblem(LinearMeasurement[] measurements) {
604 this.measurements = measurements;
605 }
606
607 public WeightedMeasurement[] getMeasurements() {
608 return measurements;
609 }
610
611 public EstimatedParameter[] getUnboundParameters() {
612 return getAllParameters();
613 }
614
615 public EstimatedParameter[] getAllParameters() {
616 HashSet<EstimatedParameter> set = new HashSet<EstimatedParameter>();
617 for (int i = 0; i < measurements.length; ++i) {
618 EstimatedParameter[] parameters = measurements[i].getParameters();
619 for (int j = 0; j < parameters.length; ++j) {
620 set.add(parameters[j]);
621 }
622 }
623 return set.toArray(new EstimatedParameter[set.size()]);
624 }
625
626 private LinearMeasurement[] measurements;
627
628 }
629
630 private static class LinearMeasurement extends WeightedMeasurement {
631
632 public LinearMeasurement(double[] factors, EstimatedParameter[] parameters,
633 double setPoint) {
634 super(1.0, setPoint);
635 this.factors = factors;
636 this.parameters = parameters;
637 }
638
639 @Override
640 public double getTheoreticalValue() {
641 double v = 0;
642 for (int i = 0; i < factors.length; ++i) {
643 v += factors[i] * parameters[i].getEstimate();
644 }
645 return v;
646 }
647
648 @Override
649 public double getPartial(EstimatedParameter parameter) {
650 for (int i = 0; i < parameters.length; ++i) {
651 if (parameters[i] == parameter) {
652 return factors[i];
653 }
654 }
655 return 0;
656 }
657
658 public EstimatedParameter[] getParameters() {
659 return parameters;
660 }
661
662 private double[] factors;
663 private EstimatedParameter[] parameters;
664 private static final long serialVersionUID = -3922448707008868580L;
665
666 }
667
668 private static class Circle implements EstimationProblem {
669
670 public Circle(double cx, double cy) {
671 this.cx = new EstimatedParameter("cx", cx);
672 this.cy = new EstimatedParameter("cy", cy);
673 points = new ArrayList<PointModel>();
674 }
675
676 public void addPoint(double px, double py) {
677 points.add(new PointModel(this, px, py));
678 }
679
680 public int getM() {
681 return points.size();
682 }
683
684 public WeightedMeasurement[] getMeasurements() {
685 return points.toArray(new PointModel[points.size()]);
686 }
687
688 public EstimatedParameter[] getAllParameters() {
689 return new EstimatedParameter[] { cx, cy };
690 }
691
692 public EstimatedParameter[] getUnboundParameters() {
693 return new EstimatedParameter[] { cx, cy };
694 }
695
696 public double getPartialRadiusX() {
697 double dRdX = 0;
698 for (PointModel point : points) {
699 dRdX += point.getPartialDiX();
700 }
701 return dRdX / points.size();
702 }
703
704 public double getPartialRadiusY() {
705 double dRdY = 0;
706 for (PointModel point : points) {
707 dRdY += point.getPartialDiY();
708 }
709 return dRdY / points.size();
710 }
711
712 public double getRadius() {
713 double r = 0;
714 for (PointModel point : points) {
715 r += point.getCenterDistance();
716 }
717 return r / points.size();
718 }
719
720 public double getX() {
721 return cx.getEstimate();
722 }
723
724 public double getY() {
725 return cy.getEstimate();
726 }
727
728 private static class PointModel extends WeightedMeasurement {
729
730 public PointModel(Circle circle, double px, double py) {
731 super(1.0, 0.0);
732 this.px = px;
733 this.py = py;
734 this.circle = circle;
735 }
736
737 @Override
738 public double getPartial(EstimatedParameter parameter) {
739 if (parameter == circle.cx) {
740 return getPartialDiX() - circle.getPartialRadiusX();
741 } else if (parameter == circle.cy) {
742 return getPartialDiY() - circle.getPartialRadiusY();
743 }
744 return 0;
745 }
746
747 public double getCenterDistance() {
748 double dx = px - circle.cx.getEstimate();
749 double dy = py - circle.cy.getEstimate();
750 return Math.sqrt(dx * dx + dy * dy);
751 }
752
753 public double getPartialDiX() {
754 return (circle.cx.getEstimate() - px) / getCenterDistance();
755 }
756
757 public double getPartialDiY() {
758 return (circle.cy.getEstimate() - py) / getCenterDistance();
759 }
760
761 @Override
762 public double getTheoreticalValue() {
763 return getCenterDistance() - circle.getRadius();
764 }
765
766 private double px;
767 private double py;
768 private transient final Circle circle;
769 private static final long serialVersionUID = 1L;
770
771 }
772
773 private EstimatedParameter cx;
774 private EstimatedParameter cy;
775 private ArrayList<PointModel> points;
776
777 }
778
779 private static class QuadraticProblem extends SimpleEstimationProblem {
780
781 private EstimatedParameter a;
782 private EstimatedParameter b;
783 private EstimatedParameter c;
784
785 public QuadraticProblem() {
786 a = new EstimatedParameter("a", 0.0);
787 b = new EstimatedParameter("b", 0.0);
788 c = new EstimatedParameter("c", 0.0);
789 addParameter(a);
790 addParameter(b);
791 addParameter(c);
792 }
793
794 public void addPoint(double x, double y, double w) {
795 addMeasurement(new LocalMeasurement(this, x, y, w));
796 }
797
798 public double theoreticalValue(double x) {
799 return ( (a.getEstimate() * x + b.getEstimate() ) * x + c.getEstimate());
800 }
801
802 private double partial(double x, EstimatedParameter parameter) {
803 if (parameter == a) {
804 return x * x;
805 } else if (parameter == b) {
806 return x;
807 } else {
808 return 1.0;
809 }
810 }
811
812 private static class LocalMeasurement extends WeightedMeasurement {
813
814 private static final long serialVersionUID = 1555043155023729130L;
815 private final double x;
816 private transient final QuadraticProblem pb;
817
818
819 public LocalMeasurement(QuadraticProblem pb, double x, double y, double w) {
820 super(w, y);
821 this.x = x;
822 this.pb = pb;
823 }
824
825 @Override
826 public double getTheoreticalValue() {
827 return pb.theoreticalValue(x);
828 }
829
830 @Override
831 public double getPartial(EstimatedParameter parameter) {
832 return pb.partial(x, parameter);
833 }
834
835 }
836 }
837
838 public static Test suite() {
839 return new TestSuite(LevenbergMarquardtEstimatorTest.class);
840 }
841
842 }