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