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  
25  import junit.framework.Test;
26  import junit.framework.TestCase;
27  import junit.framework.TestSuite;
28  
29  import org.apache.commons.math.FunctionEvaluationException;
30  import org.apache.commons.math.analysis.DifferentiableMultivariateVectorialFunction;
31  import org.apache.commons.math.analysis.MultivariateMatrixFunction;
32  import org.apache.commons.math.linear.BlockRealMatrix;
33  import org.apache.commons.math.linear.RealMatrix;
34  import org.apache.commons.math.optimization.OptimizationException;
35  import org.apache.commons.math.optimization.SimpleVectorialPointChecker;
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 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 }