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 LevenbergMarquardtEstimatorTest
091      extends TestCase {
092    
093      public LevenbergMarquardtEstimatorTest(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        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    }