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.optimization.general;
19  
20  import java.awt.geom.Point2D;
21  import java.io.Serializable;
22  import java.util.ArrayList;
23  import java.util.Arrays;
24  import java.util.List;
25  
26  import junit.framework.Test;
27  import junit.framework.TestCase;
28  import junit.framework.TestSuite;
29  
30  import org.apache.commons.math.FunctionEvaluationException;
31  import org.apache.commons.math.analysis.DifferentiableMultivariateVectorialFunction;
32  import org.apache.commons.math.analysis.MultivariateMatrixFunction;
33  import org.apache.commons.math.linear.BlockRealMatrix;
34  import org.apache.commons.math.linear.RealMatrix;
35  import org.apache.commons.math.optimization.OptimizationException;
36  import org.apache.commons.math.optimization.SimpleVectorialValueChecker;
37  import org.apache.commons.math.optimization.VectorialPointValuePair;
38  
39  /**
40   * <p>Some of the unit tests are re-implementations of the MINPACK <a
41   * href="http://www.netlib.org/minpack/ex/file17">file17</a> and <a
42   * href="http://www.netlib.org/minpack/ex/file22">file22</a> test files. 
43   * The redistribution policy for MINPACK is available <a
44   * href="http://www.netlib.org/minpack/disclaimer">here</a>, for
45   * convenience, it is reproduced below.</p>
46  
47   * <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0">
48   * <tr><td>
49   *    Minpack Copyright Notice (1999) University of Chicago.
50   *    All rights reserved
51   * </td></tr>
52   * <tr><td>
53   * Redistribution and use in source and binary forms, with or without
54   * modification, are permitted provided that the following conditions
55   * are met:
56   * <ol>
57   *  <li>Redistributions of source code must retain the above copyright
58   *      notice, this list of conditions and the following disclaimer.</li>
59   * <li>Redistributions in binary form must reproduce the above
60   *     copyright notice, this list of conditions and the following
61   *     disclaimer in the documentation and/or other materials provided
62   *     with the distribution.</li>
63   * <li>The end-user documentation included with the redistribution, if any,
64   *     must include the following acknowledgment:
65   *     <code>This product includes software developed by the University of
66   *           Chicago, as Operator of Argonne National Laboratory.</code>
67   *     Alternately, this acknowledgment may appear in the software itself,
68   *     if and wherever such third-party acknowledgments normally appear.</li>
69   * <li><strong>WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS"
70   *     WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE
71   *     UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND
72   *     THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR
73   *     IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES
74   *     OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE
75   *     OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY
76   *     OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR
77   *     USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF
78   *     THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4)
79   *     DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION
80   *     UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL
81   *     BE CORRECTED.</strong></li>
82   * <li><strong>LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT
83   *     HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF
84   *     ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT,
85   *     INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF
86   *     ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF
87   *     PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER
88   *     SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT
89   *     (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE,
90   *     EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE
91   *     POSSIBILITY OF SUCH LOSS OR DAMAGES.</strong></li>
92   * <ol></td></tr>
93   * </table>
94  
95   * @author Argonne National Laboratory. MINPACK project. March 1980 (original fortran minpack tests)
96   * @author Burton S. Garbow (original fortran minpack tests)
97   * @author Kenneth E. Hillstrom (original fortran minpack tests)
98   * @author Jorge J. More (original fortran minpack tests)
99   * @author Luc Maisonobe (non-minpack tests and minpack tests Java translation)
100  */
101 public class LevenbergMarquardtOptimizerTest
102   extends TestCase {
103 
104     public LevenbergMarquardtOptimizerTest(String name) {
105         super(name);
106     }
107 
108     public void testTrivial() throws FunctionEvaluationException, OptimizationException {
109         LinearProblem problem =
110             new LinearProblem(new double[][] { { 2 } }, new double[] { 3 });
111         LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
112         VectorialPointValuePair optimum =
113             optimizer.optimize(problem, problem.target, new double[] { 1 }, new double[] { 0 });
114         assertEquals(0, optimizer.getRMS(), 1.0e-10);
115         try {
116             optimizer.guessParametersErrors();
117             fail("an exception should have been thrown");
118         } catch (OptimizationException ee) {
119             // expected behavior
120         } catch (Exception e) {
121             fail("wrong exception caught");
122         }
123         assertEquals(1.5, optimum.getPoint()[0], 1.0e-10);
124         assertEquals(3.0, optimum.getValue()[0], 1.0e-10);
125     }
126 
127     public void testQRColumnsPermutation() throws FunctionEvaluationException, OptimizationException {
128 
129         LinearProblem problem =
130             new LinearProblem(new double[][] { { 1.0, -1.0 }, { 0.0, 2.0 }, { 1.0, -2.0 } },
131                               new double[] { 4.0, 6.0, 1.0 });
132 
133         LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
134         VectorialPointValuePair optimum =
135             optimizer.optimize(problem, problem.target, new double[] { 1, 1, 1 }, new double[] { 0, 0 });
136         assertEquals(0, optimizer.getRMS(), 1.0e-10);
137         assertEquals(7.0, optimum.getPoint()[0], 1.0e-10);
138         assertEquals(3.0, optimum.getPoint()[1], 1.0e-10);
139         assertEquals(4.0, optimum.getValue()[0], 1.0e-10);
140         assertEquals(6.0, optimum.getValue()[1], 1.0e-10);
141         assertEquals(1.0, optimum.getValue()[2], 1.0e-10);
142 
143     }
144 
145     public void testNoDependency() throws FunctionEvaluationException, OptimizationException {
146         LinearProblem problem = new LinearProblem(new double[][] {
147                 { 2, 0, 0, 0, 0, 0 },
148                 { 0, 2, 0, 0, 0, 0 },
149                 { 0, 0, 2, 0, 0, 0 },
150                 { 0, 0, 0, 2, 0, 0 },
151                 { 0, 0, 0, 0, 2, 0 },
152                 { 0, 0, 0, 0, 0, 2 }
153         }, new double[] { 0.0, 1.1, 2.2, 3.3, 4.4, 5.5 });
154         LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
155         VectorialPointValuePair optimum =
156             optimizer.optimize(problem, problem.target, new double[] { 1, 1, 1, 1, 1, 1 },
157                                new double[] { 0, 0, 0, 0, 0, 0 });
158         assertEquals(0, optimizer.getRMS(), 1.0e-10);
159         for (int i = 0; i < problem.target.length; ++i) {
160             assertEquals(0.55 * i, optimum.getPoint()[i], 1.0e-10);
161         }
162     }
163 
164     public void testOneSet() throws FunctionEvaluationException, OptimizationException {
165 
166         LinearProblem problem = new LinearProblem(new double[][] {
167                 {  1,  0, 0 },
168                 { -1,  1, 0 },
169                 {  0, -1, 1 }
170         }, new double[] { 1, 1, 1});
171         LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
172         VectorialPointValuePair optimum =
173             optimizer.optimize(problem, problem.target, new double[] { 1, 1, 1 }, new double[] { 0, 0, 0 });
174         assertEquals(0, optimizer.getRMS(), 1.0e-10);
175         assertEquals(1.0, optimum.getPoint()[0], 1.0e-10);
176         assertEquals(2.0, optimum.getPoint()[1], 1.0e-10);
177         assertEquals(3.0, optimum.getPoint()[2], 1.0e-10);
178 
179     }
180 
181     public void testTwoSets() throws FunctionEvaluationException, OptimizationException {
182         double epsilon = 1.0e-7;
183         LinearProblem problem = new LinearProblem(new double[][] {
184                 {  2,  1,   0,  4,       0, 0 },
185                 { -4, -2,   3, -7,       0, 0 },
186                 {  4,  1,  -2,  8,       0, 0 },
187                 {  0, -3, -12, -1,       0, 0 },
188                 {  0,  0,   0,  0, epsilon, 1 },
189                 {  0,  0,   0,  0,       1, 1 }
190         }, new double[] { 2, -9, 2, 2, 1 + epsilon * epsilon, 2});
191 
192         LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
193         VectorialPointValuePair optimum =
194             optimizer.optimize(problem, problem.target, new double[] { 1, 1, 1, 1, 1, 1 },
195                                new double[] { 0, 0, 0, 0, 0, 0 });
196         assertEquals(0, optimizer.getRMS(), 1.0e-10);
197         assertEquals( 3.0, optimum.getPoint()[0], 1.0e-10);
198         assertEquals( 4.0, optimum.getPoint()[1], 1.0e-10);
199         assertEquals(-1.0, optimum.getPoint()[2], 1.0e-10);
200         assertEquals(-2.0, optimum.getPoint()[3], 1.0e-10);
201         assertEquals( 1.0 + epsilon, optimum.getPoint()[4], 1.0e-10);
202         assertEquals( 1.0 - epsilon, optimum.getPoint()[5], 1.0e-10);
203 
204     }
205 
206     public void testNonInversible() throws FunctionEvaluationException, OptimizationException {
207 
208         LinearProblem problem = new LinearProblem(new double[][] {
209                 {  1, 2, -3 },
210                 {  2, 1,  3 },
211                 { -3, 0, -9 }
212         }, new double[] { 1, 1, 1 });
213  
214         LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
215         optimizer.optimize(problem, problem.target, new double[] { 1, 1, 1 }, new double[] { 0, 0, 0 });
216         assertTrue(Math.sqrt(problem.target.length) * optimizer.getRMS() > 0.6);
217         try {
218             optimizer.getCovariances();
219             fail("an exception should have been thrown");
220         } catch (OptimizationException ee) {
221             // expected behavior
222         } catch (Exception e) {
223             fail("wrong exception caught");
224         }
225 
226     }
227 
228     public void testIllConditioned() throws FunctionEvaluationException, OptimizationException {
229         LinearProblem problem1 = new LinearProblem(new double[][] {
230                 { 10.0, 7.0,  8.0,  7.0 },
231                 {  7.0, 5.0,  6.0,  5.0 },
232                 {  8.0, 6.0, 10.0,  9.0 },
233                 {  7.0, 5.0,  9.0, 10.0 }
234         }, new double[] { 32, 23, 33, 31 });
235         LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
236         VectorialPointValuePair optimum1 =
237             optimizer.optimize(problem1, problem1.target, new double[] { 1, 1, 1, 1 },
238                                new double[] { 0, 1, 2, 3 });
239         assertEquals(0, optimizer.getRMS(), 1.0e-10);
240         assertEquals(1.0, optimum1.getPoint()[0], 1.0e-10);
241         assertEquals(1.0, optimum1.getPoint()[1], 1.0e-10);
242         assertEquals(1.0, optimum1.getPoint()[2], 1.0e-10);
243         assertEquals(1.0, optimum1.getPoint()[3], 1.0e-10);
244 
245         LinearProblem problem2 = new LinearProblem(new double[][] {
246                 { 10.00, 7.00, 8.10, 7.20 },
247                 {  7.08, 5.04, 6.00, 5.00 },
248                 {  8.00, 5.98, 9.89, 9.00 },
249                 {  6.99, 4.99, 9.00, 9.98 }
250         }, new double[] { 32, 23, 33, 31 });
251         VectorialPointValuePair optimum2 =
252             optimizer.optimize(problem2, problem2.target, new double[] { 1, 1, 1, 1 },
253                                new double[] { 0, 1, 2, 3 });
254         assertEquals(0, optimizer.getRMS(), 1.0e-10);
255         assertEquals(-81.0, optimum2.getPoint()[0], 1.0e-8);
256         assertEquals(137.0, optimum2.getPoint()[1], 1.0e-8);
257         assertEquals(-34.0, optimum2.getPoint()[2], 1.0e-8);
258         assertEquals( 22.0, optimum2.getPoint()[3], 1.0e-8);
259 
260     }
261 
262     public void testMoreEstimatedParametersSimple() throws FunctionEvaluationException, OptimizationException {
263 
264         LinearProblem problem = new LinearProblem(new double[][] {
265                 { 3.0, 2.0,  0.0, 0.0 },
266                 { 0.0, 1.0, -1.0, 1.0 },
267                 { 2.0, 0.0,  1.0, 0.0 }
268         }, new double[] { 7.0, 3.0, 5.0 });
269 
270         LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
271         optimizer.optimize(problem, problem.target, new double[] { 1, 1, 1 },
272                 new double[] { 7, 6, 5, 4 });
273         assertEquals(0, optimizer.getRMS(), 1.0e-10);
274 
275     }
276 
277     public void testMoreEstimatedParametersUnsorted() throws FunctionEvaluationException, OptimizationException {
278         LinearProblem problem = new LinearProblem(new double[][] {
279                 { 1.0, 1.0,  0.0,  0.0, 0.0,  0.0 },
280                 { 0.0, 0.0,  1.0,  1.0, 1.0,  0.0 },
281                 { 0.0, 0.0,  0.0,  0.0, 1.0, -1.0 },
282                 { 0.0, 0.0, -1.0,  1.0, 0.0,  1.0 },
283                 { 0.0, 0.0,  0.0, -1.0, 1.0,  0.0 }
284        }, new double[] { 3.0, 12.0, -1.0, 7.0, 1.0 });
285 
286         LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
287         VectorialPointValuePair optimum =
288             optimizer.optimize(problem, problem.target, new double[] { 1, 1, 1, 1, 1 },
289                                new double[] { 2, 2, 2, 2, 2, 2 });
290         assertEquals(0, optimizer.getRMS(), 1.0e-10);
291         assertEquals(3.0, optimum.getPointRef()[2], 1.0e-10);
292         assertEquals(4.0, optimum.getPointRef()[3], 1.0e-10);
293         assertEquals(5.0, optimum.getPointRef()[4], 1.0e-10);
294         assertEquals(6.0, optimum.getPointRef()[5], 1.0e-10);
295 
296     }
297 
298     public void testRedundantEquations() throws FunctionEvaluationException, OptimizationException {
299         LinearProblem problem = new LinearProblem(new double[][] {
300                 { 1.0,  1.0 },
301                 { 1.0, -1.0 },
302                 { 1.0,  3.0 }
303         }, new double[] { 3.0, 1.0, 5.0 });
304 
305         LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
306         VectorialPointValuePair optimum =
307             optimizer.optimize(problem, problem.target, new double[] { 1, 1, 1 },
308                                new double[] { 1, 1 });
309         assertEquals(0, optimizer.getRMS(), 1.0e-10);
310         assertEquals(2.0, optimum.getPointRef()[0], 1.0e-10);
311         assertEquals(1.0, optimum.getPointRef()[1], 1.0e-10);
312 
313     }
314 
315     public void testInconsistentEquations() throws FunctionEvaluationException, OptimizationException {
316         LinearProblem problem = new LinearProblem(new double[][] {
317                 { 1.0,  1.0 },
318                 { 1.0, -1.0 },
319                 { 1.0,  3.0 }
320         }, new double[] { 3.0, 1.0, 4.0 });
321 
322         LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
323         optimizer.optimize(problem, problem.target, new double[] { 1, 1, 1 }, new double[] { 1, 1 });
324         assertTrue(optimizer.getRMS() > 0.1);
325 
326     }
327 
328     public void testInconsistentSizes() throws FunctionEvaluationException, OptimizationException {
329         LinearProblem problem =
330             new LinearProblem(new double[][] { { 1, 0 }, { 0, 1 } }, new double[] { -1, 1 });
331         LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
332 
333         VectorialPointValuePair optimum =
334             optimizer.optimize(problem, problem.target, new double[] { 1, 1 }, new double[] { 0, 0 });
335         assertEquals(0, optimizer.getRMS(), 1.0e-10);
336         assertEquals(-1, optimum.getPoint()[0], 1.0e-10);
337         assertEquals(+1, optimum.getPoint()[1], 1.0e-10);
338 
339         try {
340             optimizer.optimize(problem, problem.target,
341                                new double[] { 1 },
342                                new double[] { 0, 0 });
343             fail("an exception should have been thrown");
344         } catch (OptimizationException oe) {
345             // expected behavior
346         } catch (Exception e) {
347             fail("wrong exception caught");
348         }
349 
350         try {
351             optimizer.optimize(problem, new double[] { 1 },
352                                new double[] { 1 },
353                                new double[] { 0, 0 });
354             fail("an exception should have been thrown");
355         } catch (FunctionEvaluationException oe) {
356             // expected behavior
357         } catch (Exception e) {
358             fail("wrong exception caught");
359         }
360 
361     }
362 
363     public void testControlParameters() {
364         Circle circle = new Circle();
365         circle.addPoint( 30.0,  68.0);
366         circle.addPoint( 50.0,  -6.0);
367         circle.addPoint(110.0, -20.0);
368         circle.addPoint( 35.0,  15.0);
369         circle.addPoint( 45.0,  97.0);
370         checkEstimate(circle, 0.1, 10, 1.0e-14, 1.0e-16, 1.0e-10, false);
371         checkEstimate(circle, 0.1, 10, 1.0e-15, 1.0e-17, 1.0e-10, true);
372         checkEstimate(circle, 0.1,  5, 1.0e-15, 1.0e-16, 1.0e-10, true);
373         circle.addPoint(300, -300);
374         checkEstimate(circle, 0.1, 20, 1.0e-18, 1.0e-16, 1.0e-10, true);
375     }
376 
377     private void checkEstimate(DifferentiableMultivariateVectorialFunction problem,
378                                double initialStepBoundFactor, int maxCostEval,
379                                double costRelativeTolerance, double parRelativeTolerance,
380                                double orthoTolerance, boolean shouldFail) {
381         try {
382             LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
383             optimizer.setInitialStepBoundFactor(initialStepBoundFactor);
384             optimizer.setMaxIterations(maxCostEval);
385             optimizer.setCostRelativeTolerance(costRelativeTolerance);
386             optimizer.setParRelativeTolerance(parRelativeTolerance);
387             optimizer.setOrthoTolerance(orthoTolerance);
388             optimizer.optimize(problem, new double[] { 0, 0, 0, 0, 0 }, new double[] { 1, 1, 1, 1, 1 },
389                                new double[] { 98.680, 47.345 });
390             assertTrue(! shouldFail);
391         } catch (OptimizationException ee) {
392             assertTrue(shouldFail);
393         } catch (FunctionEvaluationException ee) {
394             assertTrue(shouldFail);
395         } catch (Exception e) {
396             fail("wrong exception type caught");
397         }
398     }
399 
400     public void testCircleFitting() throws FunctionEvaluationException, OptimizationException {
401         Circle circle = new Circle();
402         circle.addPoint( 30.0,  68.0);
403         circle.addPoint( 50.0,  -6.0);
404         circle.addPoint(110.0, -20.0);
405         circle.addPoint( 35.0,  15.0);
406         circle.addPoint( 45.0,  97.0);
407         LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
408         VectorialPointValuePair optimum =
409             optimizer.optimize(circle, new double[] { 0, 0, 0, 0, 0 }, new double[] { 1, 1, 1, 1, 1 },
410                                new double[] { 98.680, 47.345 });
411         assertTrue(optimizer.getEvaluations() < 10);
412         assertTrue(optimizer.getJacobianEvaluations() < 10);
413         double rms = optimizer.getRMS();
414         assertEquals(1.768262623567235,  Math.sqrt(circle.getN()) * rms,  1.0e-10);
415         Point2D.Double center = new Point2D.Double(optimum.getPointRef()[0], optimum.getPointRef()[1]);
416         assertEquals(69.96016176931406, circle.getRadius(center), 1.0e-10);
417         assertEquals(96.07590211815305, center.x,      1.0e-10);
418         assertEquals(48.13516790438953, center.y,      1.0e-10);
419         double[][] cov = optimizer.getCovariances();
420         assertEquals(1.839, cov[0][0], 0.001);
421         assertEquals(0.731, cov[0][1], 0.001);
422         assertEquals(cov[0][1], cov[1][0], 1.0e-14);
423         assertEquals(0.786, cov[1][1], 0.001);
424         double[] errors = optimizer.guessParametersErrors();
425         assertEquals(1.384, errors[0], 0.001);
426         assertEquals(0.905, errors[1], 0.001);
427 
428         // add perfect measurements and check errors are reduced
429         double  r = circle.getRadius(center);
430         for (double d= 0; d < 2 * Math.PI; d += 0.01) {
431             circle.addPoint(center.x + r * Math.cos(d), center.y + r * Math.sin(d));
432         }
433         double[] target = new double[circle.getN()];
434         Arrays.fill(target, 0.0);
435         double[] weights = new double[circle.getN()];
436         Arrays.fill(weights, 2.0);
437         optimizer.optimize(circle, target, weights, new double[] { 98.680, 47.345 });
438         cov = optimizer.getCovariances();
439         assertEquals(0.0016, cov[0][0], 0.001);
440         assertEquals(3.2e-7, cov[0][1], 1.0e-9);
441         assertEquals(cov[0][1], cov[1][0], 1.0e-14);
442         assertEquals(0.0016, cov[1][1], 0.001);
443         errors = optimizer.guessParametersErrors();
444         assertEquals(0.002, errors[0], 0.001);
445         assertEquals(0.002, errors[1], 0.001);
446 
447     }
448 
449     public void testCircleFittingBadInit() throws FunctionEvaluationException, OptimizationException {
450         Circle circle = new Circle();
451         double[][] points = new double[][] {
452                 {-0.312967,  0.072366}, {-0.339248,  0.132965}, {-0.379780,  0.202724},
453                 {-0.390426,  0.260487}, {-0.361212,  0.328325}, {-0.346039,  0.392619},
454                 {-0.280579,  0.444306}, {-0.216035,  0.470009}, {-0.149127,  0.493832},
455                 {-0.075133,  0.483271}, {-0.007759,  0.452680}, { 0.060071,  0.410235},
456                 { 0.103037,  0.341076}, { 0.118438,  0.273884}, { 0.131293,  0.192201},
457                 { 0.115869,  0.129797}, { 0.072223,  0.058396}, { 0.022884,  0.000718},
458                 {-0.053355, -0.020405}, {-0.123584, -0.032451}, {-0.216248, -0.032862},
459                 {-0.278592, -0.005008}, {-0.337655,  0.056658}, {-0.385899,  0.112526},
460                 {-0.405517,  0.186957}, {-0.415374,  0.262071}, {-0.387482,  0.343398},
461                 {-0.347322,  0.397943}, {-0.287623,  0.458425}, {-0.223502,  0.475513},
462                 {-0.135352,  0.478186}, {-0.061221,  0.483371}, { 0.003711,  0.422737},
463                 { 0.065054,  0.375830}, { 0.108108,  0.297099}, { 0.123882,  0.222850},
464                 { 0.117729,  0.134382}, { 0.085195,  0.056820}, { 0.029800, -0.019138},
465                 {-0.027520, -0.072374}, {-0.102268, -0.091555}, {-0.200299, -0.106578},
466                 {-0.292731, -0.091473}, {-0.356288, -0.051108}, {-0.420561,  0.014926},
467                 {-0.471036,  0.074716}, {-0.488638,  0.182508}, {-0.485990,  0.254068},
468                 {-0.463943,  0.338438}, {-0.406453,  0.404704}, {-0.334287,  0.466119},
469                 {-0.254244,  0.503188}, {-0.161548,  0.495769}, {-0.075733,  0.495560},
470                 { 0.001375,  0.434937}, { 0.082787,  0.385806}, { 0.115490,  0.323807},
471                 { 0.141089,  0.223450}, { 0.138693,  0.131703}, { 0.126415,  0.049174},
472                 { 0.066518, -0.010217}, {-0.005184, -0.070647}, {-0.080985, -0.103635},
473                 {-0.177377, -0.116887}, {-0.260628, -0.100258}, {-0.335756, -0.056251},
474                 {-0.405195, -0.000895}, {-0.444937,  0.085456}, {-0.484357,  0.175597},
475                 {-0.472453,  0.248681}, {-0.438580,  0.347463}, {-0.402304,  0.422428},
476                 {-0.326777,  0.479438}, {-0.247797,  0.505581}, {-0.152676,  0.519380},
477                 {-0.071754,  0.516264}, { 0.015942,  0.472802}, { 0.076608,  0.419077},
478                 { 0.127673,  0.330264}, { 0.159951,  0.262150}, { 0.153530,  0.172681},
479                 { 0.140653,  0.089229}, { 0.078666,  0.024981}, { 0.023807, -0.037022},
480                 {-0.048837, -0.077056}, {-0.127729, -0.075338}, {-0.221271, -0.067526}
481         };
482         double[] target = new double[points.length];
483         Arrays.fill(target, 0.0);
484         double[] weights = new double[points.length];
485         Arrays.fill(weights, 2.0);
486         for (int i = 0; i < points.length; ++i) {
487             circle.addPoint(points[i][0], points[i][1]);
488         }
489         LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
490         optimizer.setConvergenceChecker(new SimpleVectorialValueChecker(1.0e-10, 1.0e-10));
491         VectorialPointValuePair optimum =
492             optimizer.optimize(circle, target, weights, new double[] { -12, -12 });
493         Point2D.Double center = new Point2D.Double(optimum.getPointRef()[0], optimum.getPointRef()[1]);
494         assertTrue(optimizer.getEvaluations() < 25);
495         assertTrue(optimizer.getJacobianEvaluations() < 20);
496         assertEquals( 0.043, optimizer.getRMS(), 1.0e-3);
497         assertEquals( 0.292235,  circle.getRadius(center), 1.0e-6);
498         assertEquals(-0.151738,  center.x,      1.0e-6);
499         assertEquals( 0.2075001, center.y,      1.0e-6);
500     }
501 
502     public void testMath199() throws FunctionEvaluationException {
503         try {
504             QuadraticProblem problem = new QuadraticProblem();
505             problem.addPoint (0, -3.182591015485607);
506             problem.addPoint (1, -2.5581184967730577);
507             problem.addPoint (2, -2.1488478161387325);
508             problem.addPoint (3, -1.9122489313410047);
509             problem.addPoint (4, 1.7785661310051026);
510             new LevenbergMarquardtOptimizer().optimize(problem,
511                                                        new double[] { 0, 0, 0, 0, 0 },
512                                                        new double[] { 0.0, 4.4e-323, 1.0, 4.4e-323, 0.0 },
513                                                        new double[] { 0, 0, 0 });
514             fail("an exception should have been thrown");
515         } catch (OptimizationException ee) {
516             // expected behavior
517         }
518 
519     }
520 
521     private static class LinearProblem implements DifferentiableMultivariateVectorialFunction, Serializable {
522 
523         private static final long serialVersionUID = 703247177355019415L;
524         final RealMatrix factors;
525         final double[] target;
526         public LinearProblem(double[][] factors, double[] target) {
527             this.factors = new BlockRealMatrix(factors);
528             this.target  = target;
529         }
530 
531         public double[] value(double[] variables) {
532             return factors.operate(variables);
533         }
534 
535         public MultivariateMatrixFunction jacobian() {
536             return new MultivariateMatrixFunction() {
537                 private static final long serialVersionUID = 556396458721526234L;
538                 public double[][] value(double[] point) {
539                     return factors.getData();
540                 }
541             };
542         }
543 
544     }
545 
546     private static class Circle implements DifferentiableMultivariateVectorialFunction, Serializable {
547 
548         private static final long serialVersionUID = -4711170319243817874L;
549 
550         private ArrayList<Point2D.Double> points;
551 
552         public Circle() {
553             points  = new ArrayList<Point2D.Double>();
554         }
555 
556         public void addPoint(double px, double py) {
557             points.add(new Point2D.Double(px, py));
558         }
559 
560         public int getN() {
561             return points.size();
562         }
563 
564         public double getRadius(Point2D.Double center) {
565             double r = 0;
566             for (Point2D.Double point : points) {
567                 r += point.distance(center);
568             }
569             return r / points.size();
570         }
571 
572         private double[][] jacobian(double[] point) {
573 
574             int n = points.size();
575             Point2D.Double center = new Point2D.Double(point[0], point[1]);
576 
577             // gradient of the optimal radius
578             double dRdX = 0;
579             double dRdY = 0;
580             for (Point2D.Double pk : points) {
581                 double dk = pk.distance(center);
582                 dRdX += (center.x - pk.x) / dk;
583                 dRdY += (center.y - pk.y) / dk;
584             }
585             dRdX /= n;
586             dRdY /= n;
587 
588             // jacobian of the radius residuals
589             double[][] jacobian = new double[n][2];
590             for (int i = 0; i < n; ++i) {
591                 Point2D.Double pi = points.get(i);
592                 double di   = pi.distance(center);
593                 jacobian[i][0] = (center.x - pi.x) / di - dRdX;    
594                 jacobian[i][1] = (center.y - pi.y) / di - dRdY;    
595             }
596 
597             return jacobian;
598 
599         }
600 
601         public double[] value(double[] variables)
602         throws FunctionEvaluationException, IllegalArgumentException {
603 
604             Point2D.Double center = new Point2D.Double(variables[0], variables[1]);
605             double radius = getRadius(center);
606 
607             double[] residuals = new double[points.size()];
608             for (int i = 0; i < residuals.length; ++i) {
609                 residuals[i] = points.get(i).distance(center) - radius;
610             }
611 
612             return residuals;
613 
614         }
615 
616         public MultivariateMatrixFunction jacobian() {
617             return new MultivariateMatrixFunction() {
618                 private static final long serialVersionUID = -4340046230875165095L;
619                 public double[][] value(double[] point) {
620                     return jacobian(point);
621                 }
622             };
623         }
624 
625     }
626 
627     private static class QuadraticProblem implements DifferentiableMultivariateVectorialFunction, Serializable {
628 
629         private static final long serialVersionUID = 7072187082052755854L;
630         private List<Double> x;
631         private List<Double> y;
632 
633         public QuadraticProblem() {
634             x = new ArrayList<Double>();
635             y = new ArrayList<Double>();
636         }
637 
638         public void addPoint(double x, double y) {
639             this.x.add(x);
640             this.y.add(y);
641         }
642 
643         private double[][] jacobian(double[] variables) {
644             double[][] jacobian = new double[x.size()][3];
645             for (int i = 0; i < jacobian.length; ++i) {
646                 jacobian[i][0] = x.get(i) * x.get(i);
647                 jacobian[i][1] = x.get(i);
648                 jacobian[i][2] = 1.0;
649             }
650             return jacobian;
651         }
652 
653         public double[] value(double[] variables) {
654             double[] values = new double[x.size()];
655             for (int i = 0; i < values.length; ++i) {
656                 values[i] = (variables[0] * x.get(i) + variables[1]) * x.get(i) + variables[2];
657             }
658             return values;
659         }
660 
661         public MultivariateMatrixFunction jacobian() {
662             return new MultivariateMatrixFunction() {
663                 private static final long serialVersionUID = -8673650298627399464L;
664                 public double[][] value(double[] point) {
665                     return jacobian(point);
666                 }
667             };
668         }
669 
670     }
671 
672     public static Test suite() {
673         return new TestSuite(LevenbergMarquardtOptimizerTest.class);
674     }
675 
676 }