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