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.optimization.general;
019    
020    import java.awt.geom.Point2D;
021    import java.io.Serializable;
022    import java.util.ArrayList;
023    import java.util.Arrays;
024    import java.util.List;
025    
026    import junit.framework.Test;
027    import junit.framework.TestCase;
028    import junit.framework.TestSuite;
029    
030    import org.apache.commons.math.FunctionEvaluationException;
031    import org.apache.commons.math.analysis.DifferentiableMultivariateVectorialFunction;
032    import org.apache.commons.math.analysis.MultivariateMatrixFunction;
033    import org.apache.commons.math.linear.BlockRealMatrix;
034    import org.apache.commons.math.linear.RealMatrix;
035    import org.apache.commons.math.optimization.OptimizationException;
036    import org.apache.commons.math.optimization.SimpleVectorialValueChecker;
037    import org.apache.commons.math.optimization.VectorialPointValuePair;
038    
039    /**
040     * <p>Some of the unit tests are re-implementations of the MINPACK <a
041     * href="http://www.netlib.org/minpack/ex/file17">file17</a> and <a
042     * href="http://www.netlib.org/minpack/ex/file22">file22</a> test files. 
043     * The redistribution policy for MINPACK is available <a
044     * href="http://www.netlib.org/minpack/disclaimer">here</a>, for
045     * convenience, it is reproduced below.</p>
046    
047     * <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0">
048     * <tr><td>
049     *    Minpack Copyright Notice (1999) University of Chicago.
050     *    All rights reserved
051     * </td></tr>
052     * <tr><td>
053     * Redistribution and use in source and binary forms, with or without
054     * modification, are permitted provided that the following conditions
055     * are met:
056     * <ol>
057     *  <li>Redistributions of source code must retain the above copyright
058     *      notice, this list of conditions and the following disclaimer.</li>
059     * <li>Redistributions in binary form must reproduce the above
060     *     copyright notice, this list of conditions and the following
061     *     disclaimer in the documentation and/or other materials provided
062     *     with the distribution.</li>
063     * <li>The end-user documentation included with the redistribution, if any,
064     *     must include the following acknowledgment:
065     *     <code>This product includes software developed by the University of
066     *           Chicago, as Operator of Argonne National Laboratory.</code>
067     *     Alternately, this acknowledgment may appear in the software itself,
068     *     if and wherever such third-party acknowledgments normally appear.</li>
069     * <li><strong>WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS"
070     *     WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE
071     *     UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND
072     *     THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR
073     *     IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES
074     *     OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE
075     *     OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY
076     *     OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR
077     *     USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF
078     *     THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4)
079     *     DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION
080     *     UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL
081     *     BE CORRECTED.</strong></li>
082     * <li><strong>LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT
083     *     HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF
084     *     ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT,
085     *     INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF
086     *     ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF
087     *     PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER
088     *     SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT
089     *     (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE,
090     *     EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE
091     *     POSSIBILITY OF SUCH LOSS OR DAMAGES.</strong></li>
092     * <ol></td></tr>
093     * </table>
094    
095     * @author Argonne National Laboratory. MINPACK project. March 1980 (original fortran minpack tests)
096     * @author Burton S. Garbow (original fortran minpack tests)
097     * @author Kenneth E. Hillstrom (original fortran minpack tests)
098     * @author Jorge J. More (original fortran minpack tests)
099     * @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    }