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