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