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    
025    import junit.framework.Test;
026    import junit.framework.TestCase;
027    import junit.framework.TestSuite;
028    
029    import org.apache.commons.math.FunctionEvaluationException;
030    import org.apache.commons.math.analysis.DifferentiableMultivariateVectorialFunction;
031    import org.apache.commons.math.analysis.MultivariateMatrixFunction;
032    import org.apache.commons.math.linear.BlockRealMatrix;
033    import org.apache.commons.math.linear.RealMatrix;
034    import org.apache.commons.math.optimization.OptimizationException;
035    import org.apache.commons.math.optimization.SimpleVectorialPointChecker;
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 GaussNewtonOptimizerTest
102    extends TestCase {
103    
104        public GaussNewtonOptimizerTest(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            GaussNewtonOptimizer optimizer = new GaussNewtonOptimizer(true);
112            optimizer.setMaxIterations(100);
113            optimizer.setConvergenceChecker(new SimpleVectorialValueChecker(1.0e-6, 1.0e-6));
114            VectorialPointValuePair optimum =
115                optimizer.optimize(problem, problem.target, new double[] { 1 }, new double[] { 0 });
116            assertEquals(0, optimizer.getRMS(), 1.0e-10);
117            assertEquals(1.5, optimum.getPoint()[0], 1.0e-10);
118            assertEquals(3.0, optimum.getValue()[0], 1.0e-10);
119        }
120    
121        public void testColumnsPermutation() throws FunctionEvaluationException, OptimizationException {
122    
123            LinearProblem problem =
124                new LinearProblem(new double[][] { { 1.0, -1.0 }, { 0.0, 2.0 }, { 1.0, -2.0 } },
125                                  new double[] { 4.0, 6.0, 1.0 });
126    
127            GaussNewtonOptimizer optimizer = new GaussNewtonOptimizer(true);
128            optimizer.setMaxIterations(100);
129            optimizer.setConvergenceChecker(new SimpleVectorialValueChecker(1.0e-6, 1.0e-6));
130            VectorialPointValuePair optimum =
131                optimizer.optimize(problem, problem.target, new double[] { 1, 1, 1 }, new double[] { 0, 0 });
132            assertEquals(0, optimizer.getRMS(), 1.0e-10);
133            assertEquals(7.0, optimum.getPoint()[0], 1.0e-10);
134            assertEquals(3.0, optimum.getPoint()[1], 1.0e-10);
135            assertEquals(4.0, optimum.getValue()[0], 1.0e-10);
136            assertEquals(6.0, optimum.getValue()[1], 1.0e-10);
137            assertEquals(1.0, optimum.getValue()[2], 1.0e-10);
138    
139        }
140    
141        public void testNoDependency() throws FunctionEvaluationException, OptimizationException {
142            LinearProblem problem = new LinearProblem(new double[][] {
143                    { 2, 0, 0, 0, 0, 0 },
144                    { 0, 2, 0, 0, 0, 0 },
145                    { 0, 0, 2, 0, 0, 0 },
146                    { 0, 0, 0, 2, 0, 0 },
147                    { 0, 0, 0, 0, 2, 0 },
148                    { 0, 0, 0, 0, 0, 2 }
149            }, new double[] { 0.0, 1.1, 2.2, 3.3, 4.4, 5.5 });
150            GaussNewtonOptimizer optimizer = new GaussNewtonOptimizer(true);
151            optimizer.setMaxIterations(100);
152            optimizer.setConvergenceChecker(new SimpleVectorialValueChecker(1.0e-6, 1.0e-6));
153            VectorialPointValuePair optimum =
154                optimizer.optimize(problem, problem.target, new double[] { 1, 1, 1, 1, 1, 1 },
155                                   new double[] { 0, 0, 0, 0, 0, 0 });
156            assertEquals(0, optimizer.getRMS(), 1.0e-10);
157            for (int i = 0; i < problem.target.length; ++i) {
158                assertEquals(0.55 * i, optimum.getPoint()[i], 1.0e-10);
159            }
160        }
161    
162        public void testOneSet() throws FunctionEvaluationException, OptimizationException {
163    
164            LinearProblem problem = new LinearProblem(new double[][] {
165                    {  1,  0, 0 },
166                    { -1,  1, 0 },
167                    {  0, -1, 1 }
168            }, new double[] { 1, 1, 1});
169            GaussNewtonOptimizer optimizer = new GaussNewtonOptimizer(true);
170            optimizer.setMaxIterations(100);
171            optimizer.setConvergenceChecker(new SimpleVectorialValueChecker(1.0e-6, 1.0e-6));
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            GaussNewtonOptimizer optimizer = new GaussNewtonOptimizer(true);
193            optimizer.setMaxIterations(100);
194            optimizer.setConvergenceChecker(new SimpleVectorialValueChecker(1.0e-6, 1.0e-6));
195            VectorialPointValuePair optimum =
196                optimizer.optimize(problem, problem.target, new double[] { 1, 1, 1, 1, 1, 1 },
197                                   new double[] { 0, 0, 0, 0, 0, 0 });
198            assertEquals(0, optimizer.getRMS(), 1.0e-10);
199            assertEquals( 3.0, optimum.getPoint()[0], 1.0e-10);
200            assertEquals( 4.0, optimum.getPoint()[1], 1.0e-10);
201            assertEquals(-1.0, optimum.getPoint()[2], 1.0e-10);
202            assertEquals(-2.0, optimum.getPoint()[3], 1.0e-10);
203            assertEquals( 1.0 + epsilon, optimum.getPoint()[4], 1.0e-10);
204            assertEquals( 1.0 - epsilon, optimum.getPoint()[5], 1.0e-10);
205    
206        }
207    
208        public void testNonInversible() {
209    
210            LinearProblem problem = new LinearProblem(new double[][] {
211                    {  1, 2, -3 },
212                    {  2, 1,  3 },
213                    { -3, 0, -9 }
214            }, new double[] { 1, 1, 1 });
215            GaussNewtonOptimizer optimizer = new GaussNewtonOptimizer(true);
216            optimizer.setMaxIterations(100);
217            optimizer.setConvergenceChecker(new SimpleVectorialValueChecker(1.0e-6, 1.0e-6));
218            try {
219                optimizer.optimize(problem, problem.target, new double[] { 1, 1, 1 }, new double[] { 0, 0, 0 });
220                fail("an exception should have been caught");
221            } catch (OptimizationException ee) {
222                // expected behavior
223            } catch (Exception e) {
224                fail("wrong exception type caught");
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            GaussNewtonOptimizer optimizer = new GaussNewtonOptimizer(true);
236            optimizer.setMaxIterations(100);
237            optimizer.setConvergenceChecker(new SimpleVectorialValueChecker(1.0e-6, 1.0e-6));
238            VectorialPointValuePair optimum1 =
239                optimizer.optimize(problem1, problem1.target, new double[] { 1, 1, 1, 1 },
240                                   new double[] { 0, 1, 2, 3 });
241            assertEquals(0, optimizer.getRMS(), 1.0e-10);
242            assertEquals(1.0, optimum1.getPoint()[0], 1.0e-10);
243            assertEquals(1.0, optimum1.getPoint()[1], 1.0e-10);
244            assertEquals(1.0, optimum1.getPoint()[2], 1.0e-10);
245            assertEquals(1.0, optimum1.getPoint()[3], 1.0e-10);
246    
247            LinearProblem problem2 = new LinearProblem(new double[][] {
248                    { 10.00, 7.00, 8.10, 7.20 },
249                    {  7.08, 5.04, 6.00, 5.00 },
250                    {  8.00, 5.98, 9.89, 9.00 },
251                    {  6.99, 4.99, 9.00, 9.98 }
252            }, new double[] { 32, 23, 33, 31 });
253            VectorialPointValuePair optimum2 =
254                optimizer.optimize(problem2, problem2.target, new double[] { 1, 1, 1, 1 },
255                                   new double[] { 0, 1, 2, 3 });
256            assertEquals(0, optimizer.getRMS(), 1.0e-10);
257            assertEquals(-81.0, optimum2.getPoint()[0], 1.0e-8);
258            assertEquals(137.0, optimum2.getPoint()[1], 1.0e-8);
259            assertEquals(-34.0, optimum2.getPoint()[2], 1.0e-8);
260            assertEquals( 22.0, optimum2.getPoint()[3], 1.0e-8);
261    
262        }
263    
264        public void testMoreEstimatedParametersSimple() {
265    
266            LinearProblem problem = new LinearProblem(new double[][] {
267                    { 3.0, 2.0,  0.0, 0.0 },
268                    { 0.0, 1.0, -1.0, 1.0 },
269                    { 2.0, 0.0,  1.0, 0.0 }
270            }, new double[] { 7.0, 3.0, 5.0 });
271    
272            GaussNewtonOptimizer optimizer = new GaussNewtonOptimizer(true);
273            optimizer.setMaxIterations(100);
274            optimizer.setConvergenceChecker(new SimpleVectorialValueChecker(1.0e-6, 1.0e-6));
275            try {
276                optimizer.optimize(problem, problem.target, new double[] { 1, 1, 1 },
277                                   new double[] { 7, 6, 5, 4 });
278                fail("an exception should have been caught");
279            } catch (OptimizationException ee) {
280                // expected behavior
281            } catch (Exception e) {
282                fail("wrong exception type caught");
283            }
284    
285        }
286    
287        public void testMoreEstimatedParametersUnsorted() {
288            LinearProblem problem = new LinearProblem(new double[][] {
289                     { 1.0, 1.0,  0.0,  0.0, 0.0,  0.0 },
290                     { 0.0, 0.0,  1.0,  1.0, 1.0,  0.0 },
291                     { 0.0, 0.0,  0.0,  0.0, 1.0, -1.0 },
292                     { 0.0, 0.0, -1.0,  1.0, 0.0,  1.0 },
293                     { 0.0, 0.0,  0.0, -1.0, 1.0,  0.0 }
294            }, new double[] { 3.0, 12.0, -1.0, 7.0, 1.0 });
295            GaussNewtonOptimizer optimizer = new GaussNewtonOptimizer(true);
296            optimizer.setMaxIterations(100);
297            optimizer.setConvergenceChecker(new SimpleVectorialValueChecker(1.0e-6, 1.0e-6));
298            try {
299                optimizer.optimize(problem, problem.target, new double[] { 1, 1, 1, 1, 1 },
300                                   new double[] { 2, 2, 2, 2, 2, 2 });
301                fail("an exception should have been caught");
302            } catch (OptimizationException ee) {
303                // expected behavior
304            } catch (Exception e) {
305                fail("wrong exception type caught");
306            }
307        }
308    
309        public void testRedundantEquations() throws FunctionEvaluationException, OptimizationException {
310            LinearProblem problem = new LinearProblem(new double[][] {
311                    { 1.0,  1.0 },
312                    { 1.0, -1.0 },
313                    { 1.0,  3.0 }
314            }, new double[] { 3.0, 1.0, 5.0 });
315    
316            GaussNewtonOptimizer optimizer = new GaussNewtonOptimizer(true);
317            optimizer.setMaxIterations(100);
318            optimizer.setConvergenceChecker(new SimpleVectorialValueChecker(1.0e-6, 1.0e-6));
319            VectorialPointValuePair optimum =
320                optimizer.optimize(problem, problem.target, new double[] { 1, 1, 1 },
321                                   new double[] { 1, 1 });
322            assertEquals(0, optimizer.getRMS(), 1.0e-10);
323            assertEquals(2.0, optimum.getPoint()[0], 1.0e-8);
324            assertEquals(1.0, optimum.getPoint()[1], 1.0e-8);
325    
326        }
327    
328        public void testInconsistentEquations() throws FunctionEvaluationException, OptimizationException {
329            LinearProblem problem = new LinearProblem(new double[][] {
330                    { 1.0,  1.0 },
331                    { 1.0, -1.0 },
332                    { 1.0,  3.0 }
333            }, new double[] { 3.0, 1.0, 4.0 });
334    
335            GaussNewtonOptimizer optimizer = new GaussNewtonOptimizer(true);
336            optimizer.setMaxIterations(100);
337            optimizer.setConvergenceChecker(new SimpleVectorialValueChecker(1.0e-6, 1.0e-6));
338            optimizer.optimize(problem, problem.target, new double[] { 1, 1, 1 }, new double[] { 1, 1 });
339            assertTrue(optimizer.getRMS() > 0.1);
340    
341        }
342    
343        public void testInconsistentSizes() throws FunctionEvaluationException, OptimizationException {
344            LinearProblem problem =
345                new LinearProblem(new double[][] { { 1, 0 }, { 0, 1 } }, new double[] { -1, 1 });
346            GaussNewtonOptimizer optimizer = new GaussNewtonOptimizer(true);
347            optimizer.setMaxIterations(100);
348            optimizer.setConvergenceChecker(new SimpleVectorialValueChecker(1.0e-6, 1.0e-6));
349    
350            VectorialPointValuePair optimum =
351                optimizer.optimize(problem, problem.target, new double[] { 1, 1 }, new double[] { 0, 0 });
352            assertEquals(0, optimizer.getRMS(), 1.0e-10);
353            assertEquals(-1, optimum.getPoint()[0], 1.0e-10);
354            assertEquals(+1, optimum.getPoint()[1], 1.0e-10);
355    
356            try {
357                optimizer.optimize(problem, problem.target,
358                                   new double[] { 1 },
359                                   new double[] { 0, 0 });
360                fail("an exception should have been thrown");
361            } catch (OptimizationException oe) {
362                // expected behavior
363            } catch (Exception e) {
364                fail("wrong exception caught");
365            }
366    
367            try {
368                optimizer.optimize(problem, new double[] { 1 },
369                                   new double[] { 1 },
370                                   new double[] { 0, 0 });
371                fail("an exception should have been thrown");
372            } catch (FunctionEvaluationException oe) {
373                // expected behavior
374            } catch (Exception e) {
375                fail("wrong exception caught");
376            }
377    
378        }
379    
380        public void testMaxIterations() {
381            Circle circle = new Circle();
382            circle.addPoint( 30.0,  68.0);
383            circle.addPoint( 50.0,  -6.0);
384            circle.addPoint(110.0, -20.0);
385            circle.addPoint( 35.0,  15.0);
386            circle.addPoint( 45.0,  97.0);
387            GaussNewtonOptimizer optimizer = new GaussNewtonOptimizer(true);
388            optimizer.setMaxIterations(100);
389            optimizer.setConvergenceChecker(new SimpleVectorialPointChecker(1.0e-30, 1.0e-30));
390            try {
391                optimizer.optimize(circle, new double[] { 0, 0, 0, 0, 0 },
392                                   new double[] { 1, 1, 1, 1, 1 },
393                                   new double[] { 98.680, 47.345 });
394                fail("an exception should have been caught");
395            } catch (OptimizationException ee) {
396                // expected behavior
397            } catch (Exception e) {
398                fail("wrong exception type caught");
399            }
400        }
401    
402        public void testCircleFitting() throws FunctionEvaluationException, OptimizationException {
403            Circle circle = new Circle();
404            circle.addPoint( 30.0,  68.0);
405            circle.addPoint( 50.0,  -6.0);
406            circle.addPoint(110.0, -20.0);
407            circle.addPoint( 35.0,  15.0);
408            circle.addPoint( 45.0,  97.0);
409            GaussNewtonOptimizer optimizer = new GaussNewtonOptimizer(true);
410            optimizer.setMaxIterations(100);
411            optimizer.setConvergenceChecker(new SimpleVectorialValueChecker(1.0e-13, 1.0e-13));
412            VectorialPointValuePair optimum =
413                optimizer.optimize(circle, new double[] { 0, 0, 0, 0, 0 },
414                                   new double[] { 1, 1, 1, 1, 1 },
415                                   new double[] { 98.680, 47.345 });
416            assertEquals(1.768262623567235,  Math.sqrt(circle.getN()) * optimizer.getRMS(),  1.0e-10);
417            Point2D.Double center = new Point2D.Double(optimum.getPointRef()[0], optimum.getPointRef()[1]);
418            assertEquals(69.96016175359975, circle.getRadius(center), 1.0e-10);
419            assertEquals(96.07590209601095, center.x, 1.0e-10);
420            assertEquals(48.135167894714,   center.y, 1.0e-10);
421        }
422    
423        public void testCircleFittingBadInit() throws FunctionEvaluationException, OptimizationException {
424            Circle circle = new Circle();
425            double[][] points = new double[][] {
426                    {-0.312967,  0.072366}, {-0.339248,  0.132965}, {-0.379780,  0.202724},
427                    {-0.390426,  0.260487}, {-0.361212,  0.328325}, {-0.346039,  0.392619},
428                    {-0.280579,  0.444306}, {-0.216035,  0.470009}, {-0.149127,  0.493832},
429                    {-0.075133,  0.483271}, {-0.007759,  0.452680}, { 0.060071,  0.410235},
430                    { 0.103037,  0.341076}, { 0.118438,  0.273884}, { 0.131293,  0.192201},
431                    { 0.115869,  0.129797}, { 0.072223,  0.058396}, { 0.022884,  0.000718},
432                    {-0.053355, -0.020405}, {-0.123584, -0.032451}, {-0.216248, -0.032862},
433                    {-0.278592, -0.005008}, {-0.337655,  0.056658}, {-0.385899,  0.112526},
434                    {-0.405517,  0.186957}, {-0.415374,  0.262071}, {-0.387482,  0.343398},
435                    {-0.347322,  0.397943}, {-0.287623,  0.458425}, {-0.223502,  0.475513},
436                    {-0.135352,  0.478186}, {-0.061221,  0.483371}, { 0.003711,  0.422737},
437                    { 0.065054,  0.375830}, { 0.108108,  0.297099}, { 0.123882,  0.222850},
438                    { 0.117729,  0.134382}, { 0.085195,  0.056820}, { 0.029800, -0.019138},
439                    {-0.027520, -0.072374}, {-0.102268, -0.091555}, {-0.200299, -0.106578},
440                    {-0.292731, -0.091473}, {-0.356288, -0.051108}, {-0.420561,  0.014926},
441                    {-0.471036,  0.074716}, {-0.488638,  0.182508}, {-0.485990,  0.254068},
442                    {-0.463943,  0.338438}, {-0.406453,  0.404704}, {-0.334287,  0.466119},
443                    {-0.254244,  0.503188}, {-0.161548,  0.495769}, {-0.075733,  0.495560},
444                    { 0.001375,  0.434937}, { 0.082787,  0.385806}, { 0.115490,  0.323807},
445                    { 0.141089,  0.223450}, { 0.138693,  0.131703}, { 0.126415,  0.049174},
446                    { 0.066518, -0.010217}, {-0.005184, -0.070647}, {-0.080985, -0.103635},
447                    {-0.177377, -0.116887}, {-0.260628, -0.100258}, {-0.335756, -0.056251},
448                    {-0.405195, -0.000895}, {-0.444937,  0.085456}, {-0.484357,  0.175597},
449                    {-0.472453,  0.248681}, {-0.438580,  0.347463}, {-0.402304,  0.422428},
450                    {-0.326777,  0.479438}, {-0.247797,  0.505581}, {-0.152676,  0.519380},
451                    {-0.071754,  0.516264}, { 0.015942,  0.472802}, { 0.076608,  0.419077},
452                    { 0.127673,  0.330264}, { 0.159951,  0.262150}, { 0.153530,  0.172681},
453                    { 0.140653,  0.089229}, { 0.078666,  0.024981}, { 0.023807, -0.037022},
454                    {-0.048837, -0.077056}, {-0.127729, -0.075338}, {-0.221271, -0.067526}
455            };
456            double[] target = new double[points.length];
457            Arrays.fill(target, 0.0);
458            double[] weights = new double[points.length];
459            Arrays.fill(weights, 2.0);
460            for (int i = 0; i < points.length; ++i) {
461                circle.addPoint(points[i][0], points[i][1]);
462            }
463            GaussNewtonOptimizer optimizer = new GaussNewtonOptimizer(true);
464            optimizer.setMaxIterations(100);
465            optimizer.setConvergenceChecker(new SimpleVectorialValueChecker(1.0e-6, 1.0e-6));
466            try {
467                optimizer.optimize(circle, target, weights, new double[] { -12, -12 });
468                fail("an exception should have been caught");
469            } catch (OptimizationException ee) {
470                // expected behavior
471            } catch (Exception e) {
472                fail("wrong exception type caught");
473            }
474    
475            VectorialPointValuePair optimum =
476                optimizer.optimize(circle, target, weights, new double[] { 0, 0 });
477            assertEquals(-0.1517383071957963, optimum.getPointRef()[0], 1.0e-8);
478            assertEquals(0.2074999736353867,  optimum.getPointRef()[1], 1.0e-8);
479            assertEquals(0.04268731682389561, optimizer.getRMS(),       1.0e-8);
480    
481        }
482    
483        private static class LinearProblem implements DifferentiableMultivariateVectorialFunction, Serializable {
484    
485            private static final long serialVersionUID = -8804268799379350190L;
486            final RealMatrix factors;
487            final double[] target;
488            public LinearProblem(double[][] factors, double[] target) {
489                this.factors = new BlockRealMatrix(factors);
490                this.target  = target;
491            }
492    
493            public double[] value(double[] variables) {
494                return factors.operate(variables);
495            }
496    
497            public MultivariateMatrixFunction jacobian() {
498                return new MultivariateMatrixFunction() {
499                    private static final long serialVersionUID = -8387467946663627585L;
500                    public double[][] value(double[] point) {
501                        return factors.getData();
502                    }
503                };
504            }
505    
506        }
507    
508        private static class Circle implements DifferentiableMultivariateVectorialFunction, Serializable {
509    
510            private static final long serialVersionUID = -7165774454925027042L;
511            private ArrayList<Point2D.Double> points;
512    
513            public Circle() {
514                points  = new ArrayList<Point2D.Double>();
515            }
516    
517            public void addPoint(double px, double py) {
518                points.add(new Point2D.Double(px, py));
519            }
520    
521            public int getN() {
522                return points.size();
523            }
524    
525            public double getRadius(Point2D.Double center) {
526                double r = 0;
527                for (Point2D.Double point : points) {
528                    r += point.distance(center);
529                }
530                return r / points.size();
531            }
532    
533            private double[][] jacobian(double[] variables) {
534    
535                int n = points.size();
536                Point2D.Double center = new Point2D.Double(variables[0], variables[1]);
537    
538                // gradient of the optimal radius
539                double dRdX = 0;
540                double dRdY = 0;
541                for (Point2D.Double pk : points) {
542                    double dk = pk.distance(center);
543                    dRdX += (center.x - pk.x) / dk;
544                    dRdY += (center.y - pk.y) / dk;
545                }
546                dRdX /= n;
547                dRdY /= n;
548    
549                // jacobian of the radius residuals
550                double[][] jacobian = new double[n][2];
551                for (int i = 0; i < n; ++i) {
552                    Point2D.Double pi = points.get(i);
553                    double di   = pi.distance(center);
554                    jacobian[i][0] = (center.x - pi.x) / di - dRdX;    
555                    jacobian[i][1] = (center.y - pi.y) / di - dRdY;    
556               }
557    
558                return jacobian;
559    
560            }
561    
562            public double[] value(double[] variables) {
563    
564                Point2D.Double center = new Point2D.Double(variables[0], variables[1]);
565                double radius = getRadius(center);
566    
567                double[] residuals = new double[points.size()];
568                for (int i = 0; i < residuals.length; ++i) {
569                    residuals[i] = points.get(i).distance(center) - radius;
570                }
571    
572                return residuals;
573    
574            }
575    
576            public MultivariateMatrixFunction jacobian() {
577                return new MultivariateMatrixFunction() {
578                    private static final long serialVersionUID = -4340046230875165095L;
579                    public double[][] value(double[] point) {
580                        return jacobian(point);
581                    }
582                };
583            }
584    
585        }
586    
587        public static Test suite() {
588            return new TestSuite(GaussNewtonOptimizerTest.class);
589        }
590    
591    }