001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    
018    package org.apache.commons.math.estimation;
019    
020    import java.util.ArrayList;
021    import java.util.HashSet;
022    
023    import junit.framework.Test;
024    import junit.framework.TestCase;
025    import junit.framework.TestSuite;
026    
027    /**
028     * <p>Some of the unit tests are re-implementations of the MINPACK <a
029     * href="http://www.netlib.org/minpack/ex/file17">file17</a> and <a
030     * href="http://www.netlib.org/minpack/ex/file22">file22</a> test files. 
031     * The redistribution policy for MINPACK is available <a
032     * href="http://www.netlib.org/minpack/disclaimer">here</a>, for
033     * convenience, it is reproduced below.</p>
034    
035     * <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0">
036     * <tr><td>
037     *    Minpack Copyright Notice (1999) University of Chicago.
038     *    All rights reserved
039     * </td></tr>
040     * <tr><td>
041     * Redistribution and use in source and binary forms, with or without
042     * modification, are permitted provided that the following conditions
043     * are met:
044     * <ol>
045     *  <li>Redistributions of source code must retain the above copyright
046     *      notice, this list of conditions and the following disclaimer.</li>
047     * <li>Redistributions in binary form must reproduce the above
048     *     copyright notice, this list of conditions and the following
049     *     disclaimer in the documentation and/or other materials provided
050     *     with the distribution.</li>
051     * <li>The end-user documentation included with the redistribution, if any,
052     *     must include the following acknowledgment:
053     *     <code>This product includes software developed by the University of
054     *           Chicago, as Operator of Argonne National Laboratory.</code>
055     *     Alternately, this acknowledgment may appear in the software itself,
056     *     if and wherever such third-party acknowledgments normally appear.</li>
057     * <li><strong>WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS"
058     *     WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE
059     *     UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND
060     *     THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR
061     *     IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES
062     *     OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE
063     *     OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY
064     *     OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR
065     *     USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF
066     *     THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4)
067     *     DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION
068     *     UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL
069     *     BE CORRECTED.</strong></li>
070     * <li><strong>LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT
071     *     HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF
072     *     ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT,
073     *     INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF
074     *     ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF
075     *     PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER
076     *     SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT
077     *     (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE,
078     *     EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE
079     *     POSSIBILITY OF SUCH LOSS OR DAMAGES.</strong></li>
080     * <ol></td></tr>
081     * </table>
082    
083     * @author Argonne National Laboratory. MINPACK project. March 1980 (original fortran minpack tests)
084     * @author Burton S. Garbow (original fortran minpack tests)
085     * @author Kenneth E. Hillstrom (original fortran minpack tests)
086     * @author Jorge J. More (original fortran minpack tests)
087     * @author Luc Maisonobe (non-minpack tests and minpack tests Java translation)
088     */
089    @Deprecated
090    public class GaussNewtonEstimatorTest
091      extends TestCase {
092    
093      public GaussNewtonEstimatorTest(String name) {
094        super(name);
095      }
096    
097      public void testTrivial() throws EstimationException {
098        LinearProblem problem =
099          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    }