1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
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   * <p>Some of the unit tests are re-implementations of the MINPACK <a
29   * href="http://www.netlib.org/minpack/ex/file17">file17</a> and <a
30   * href="http://www.netlib.org/minpack/ex/file22">file22</a> test files. 
31   * The redistribution policy for MINPACK is available <a
32   * href="http://www.netlib.org/minpack/disclaimer">here</a>, for
33   * convenience, it is reproduced below.</p>
34  
35   * <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0">
36   * <tr><td>
37   *    Minpack Copyright Notice (1999) University of Chicago.
38   *    All rights reserved
39   * </td></tr>
40   * <tr><td>
41   * Redistribution and use in source and binary forms, with or without
42   * modification, are permitted provided that the following conditions
43   * are met:
44   * <ol>
45   *  <li>Redistributions of source code must retain the above copyright
46   *      notice, this list of conditions and the following disclaimer.</li>
47   * <li>Redistributions in binary form must reproduce the above
48   *     copyright notice, this list of conditions and the following
49   *     disclaimer in the documentation and/or other materials provided
50   *     with the distribution.</li>
51   * <li>The end-user documentation included with the redistribution, if any,
52   *     must include the following acknowledgment:
53   *     <code>This product includes software developed by the University of
54   *           Chicago, as Operator of Argonne National Laboratory.</code>
55   *     Alternately, this acknowledgment may appear in the software itself,
56   *     if and wherever such third-party acknowledgments normally appear.</li>
57   * <li><strong>WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS"
58   *     WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE
59   *     UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND
60   *     THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR
61   *     IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES
62   *     OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE
63   *     OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY
64   *     OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR
65   *     USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF
66   *     THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4)
67   *     DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION
68   *     UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL
69   *     BE CORRECTED.</strong></li>
70   * <li><strong>LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT
71   *     HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF
72   *     ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT,
73   *     INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF
74   *     ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF
75   *     PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER
76   *     SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT
77   *     (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE,
78   *     EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE
79   *     POSSIBILITY OF SUCH LOSS OR DAMAGES.</strong></li>
80   * <ol></td></tr>
81   * </table>
82  
83   * @author Argonne National Laboratory. MINPACK project. March 1980 (original fortran minpack tests)
84   * @author Burton S. Garbow (original fortran minpack tests)
85   * @author Kenneth E. Hillstrom (original fortran minpack tests)
86   * @author Jorge J. More (original fortran minpack tests)
87   * @author Luc Maisonobe (non-minpack tests and minpack tests Java translation)
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       // 4 elements sub-problem
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       // 2 elements sub-problem
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       // expected behavior
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         // expected behavior
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         // expected behavior
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         // expected behavior
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         // expected behavior
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 }