commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r754727 [2/3] - in /commons/proper/math/trunk/src: java/org/apache/commons/math/optimization/ java/org/apache/commons/math/optimization/general/ test/org/apache/commons/math/optimization/general/
Date Sun, 15 Mar 2009 19:11:06 GMT
Copied: commons/proper/math/trunk/src/test/org/apache/commons/math/optimization/general/LevenbergMarquardtOptimizerTest.java (from r753644, commons/proper/math/trunk/src/test/org/apache/commons/math/optimization/general/LevenbergMarquardtEstimatorTest.java)
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/org/apache/commons/math/optimization/general/LevenbergMarquardtOptimizerTest.java?p2=commons/proper/math/trunk/src/test/org/apache/commons/math/optimization/general/LevenbergMarquardtOptimizerTest.java&p1=commons/proper/math/trunk/src/test/org/apache/commons/math/optimization/general/LevenbergMarquardtEstimatorTest.java&r1=753644&r2=754727&rev=754727&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/org/apache/commons/math/optimization/general/LevenbergMarquardtEstimatorTest.java (original)
+++ commons/proper/math/trunk/src/test/org/apache/commons/math/optimization/general/LevenbergMarquardtOptimizerTest.java Sun Mar 15 19:11:02 2009
@@ -17,16 +17,23 @@
 
 package org.apache.commons.math.optimization.general;
 
+import java.awt.geom.Point2D;
 import java.util.ArrayList;
-import java.util.HashSet;
-
-import org.apache.commons.math.optimization.OptimizationException;
-
+import java.util.Arrays;
+import java.util.List;
 
 import junit.framework.Test;
 import junit.framework.TestCase;
 import junit.framework.TestSuite;
 
+import org.apache.commons.math.linear.DenseRealMatrix;
+import org.apache.commons.math.linear.RealMatrix;
+import org.apache.commons.math.optimization.ObjectiveException;
+import org.apache.commons.math.optimization.OptimizationException;
+import org.apache.commons.math.optimization.SimpleVectorialValueChecker;
+import org.apache.commons.math.optimization.VectorialDifferentiableObjectiveFunction;
+import org.apache.commons.math.optimization.VectorialPointValuePair;
+
 /**
  * <p>Some of the unit tests are re-implementations of the MINPACK <a
  * href="http://www.netlib.org/minpack/ex/file17">file17</a> and <a
@@ -89,758 +96,554 @@
  * @author Jorge J. More (original fortran minpack tests)
  * @author Luc Maisonobe (non-minpack tests and minpack tests Java translation)
  */
-public class LevenbergMarquardtEstimatorTest
+public class LevenbergMarquardtOptimizerTest
   extends TestCase {
 
-  public LevenbergMarquardtEstimatorTest(String name) {
-    super(name);
-  }
-
-  public void testTrivial() throws OptimizationException {
-    LinearProblem problem =
-      new LinearProblem(new LinearMeasurement[] {
-        new LinearMeasurement(new double[] {2},
-                              new EstimatedParameter[] {
-                                 new EstimatedParameter("p0", 0)
-                              }, 3.0)
-      });
-    LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
-    estimator.estimate(problem);
-    assertEquals(0, estimator.getRMS(problem), 1.0e-10);
-    try {
-        estimator.guessParametersErrors(problem);
-        fail("an exception should have been thrown");
-    } catch (OptimizationException ee) {
-        // expected behavior
-    } catch (Exception e) {
-        fail("wrong exception caught");
-    }
-    assertEquals(1.5,
-                 problem.getUnboundParameters()[0].getEstimate(),
-                 1.0e-10);
-   }
-
-  public void testQRColumnsPermutation() throws OptimizationException {
-
-    EstimatedParameter[] x = {
-       new EstimatedParameter("p0", 0), new EstimatedParameter("p1", 0)
-    };
-    LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
-      new LinearMeasurement(new double[] { 1.0, -1.0 },
-                            new EstimatedParameter[] { x[0], x[1] },
-                            4.0),
-      new LinearMeasurement(new double[] { 2.0 },
-                            new EstimatedParameter[] { x[1] },
-                            6.0),
-      new LinearMeasurement(new double[] { 1.0, -2.0 },
-                            new EstimatedParameter[] { x[0], x[1] },
-                            1.0)
-    });
-
-    LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
-    estimator.estimate(problem);
-    assertEquals(0, estimator.getRMS(problem), 1.0e-10);
-    assertEquals(7.0, x[0].getEstimate(), 1.0e-10);
-    assertEquals(3.0, x[1].getEstimate(), 1.0e-10);
-
-  }
-
-  public void testNoDependency() throws OptimizationException {
-    EstimatedParameter[] p = new EstimatedParameter[] {
-      new EstimatedParameter("p0", 0),
-      new EstimatedParameter("p1", 0),
-      new EstimatedParameter("p2", 0),
-      new EstimatedParameter("p3", 0),
-      new EstimatedParameter("p4", 0),
-      new EstimatedParameter("p5", 0)
-    };
-    LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
-      new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[0] }, 0.0),
-      new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[1] }, 1.1),
-      new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[2] }, 2.2),
-      new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[3] }, 3.3),
-      new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[4] }, 4.4),
-      new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[5] }, 5.5)
-    });
-  LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
-  estimator.estimate(problem);
-  assertEquals(0, estimator.getRMS(problem), 1.0e-10);
-  for (int i = 0; i < p.length; ++i) {
-    assertEquals(0.55 * i, p[i].getEstimate(), 1.0e-10);
-  }
-}
+    public LevenbergMarquardtOptimizerTest(String name) {
+        super(name);
+    }
+
+    public void testTrivial() throws ObjectiveException, OptimizationException {
+        LinearProblem problem =
+            new LinearProblem(new double[][] { { 2 } }, new double[] { 3 });
+        LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
+        VectorialPointValuePair optimum =
+            optimizer.optimize(problem, problem.target, new double[] { 1 }, new double[] { 0 });
+        assertEquals(0, optimizer.getRMS(), 1.0e-10);
+        try {
+            optimizer.guessParametersErrors();
+            fail("an exception should have been thrown");
+        } catch (OptimizationException ee) {
+            // expected behavior
+        } catch (Exception e) {
+            fail("wrong exception caught");
+        }
+        assertEquals(1.5, optimum.getPoint()[0], 1.0e-10);
+    }
+
+    public void testQRColumnsPermutation() throws ObjectiveException, OptimizationException {
+
+        LinearProblem problem =
+            new LinearProblem(new double[][] { { 1.0, -1.0 }, { 0.0, 2.0 }, { 1.0, -2.0 } },
+                    new double[] { 4.0, 6.0, 1.0 });
+
+        LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
+        VectorialPointValuePair optimum =
+            optimizer.optimize(problem, problem.target, new double[] { 1, 1, 1 }, new double[] { 0, 0 });
+        assertEquals(0, optimizer.getRMS(), 1.0e-10);
+        assertEquals(7.0, optimum.getPoint()[0], 1.0e-10);
+        assertEquals(3.0, optimum.getPoint()[1], 1.0e-10);
+
+    }
+
+    public void testNoDependency() throws ObjectiveException, OptimizationException {
+        LinearProblem problem = new LinearProblem(new double[][] {
+                { 2, 0, 0, 0, 0, 0 },
+                { 0, 2, 0, 0, 0, 0 },
+                { 0, 0, 2, 0, 0, 0 },
+                { 0, 0, 0, 2, 0, 0 },
+                { 0, 0, 0, 0, 2, 0 },
+                { 0, 0, 0, 0, 0, 2 }
+        }, new double[] { 0.0, 1.1, 2.2, 3.3, 4.4, 5.5 });
+        LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
+        VectorialPointValuePair optimum =
+            optimizer.optimize(problem, problem.target, new double[] { 1, 1, 1, 1, 1, 1 },
+                               new double[] { 0, 0, 0, 0, 0, 0 });
+        assertEquals(0, optimizer.getRMS(), 1.0e-10);
+        for (int i = 0; i < problem.target.length; ++i) {
+            assertEquals(0.55 * i, optimum.getPoint()[i], 1.0e-10);
+        }
+    }
+
+    public void testOneSet() throws ObjectiveException, OptimizationException {
+
+        LinearProblem problem = new LinearProblem(new double[][] {
+                {  1,  0, 0 },
+                { -1,  1, 0 },
+                {  0, -1, 1 }
+        }, new double[] { 1, 1, 1});
+        LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
+        VectorialPointValuePair optimum =
+            optimizer.optimize(problem, problem.target, new double[] { 1, 1, 1 }, new double[] { 0, 0, 0 });
+        assertEquals(0, optimizer.getRMS(), 1.0e-10);
+        assertEquals(1.0, optimum.getPoint()[0], 1.0e-10);
+        assertEquals(2.0, optimum.getPoint()[1], 1.0e-10);
+        assertEquals(3.0, optimum.getPoint()[2], 1.0e-10);
+
+    }
+
+    public void testTwoSets() throws ObjectiveException, OptimizationException {
+        double epsilon = 1.0e-7;
+        LinearProblem problem = new LinearProblem(new double[][] {
+                {  2,  1,   0,  4,       0, 0 },
+                { -4, -2,   3, -7,       0, 0 },
+                {  4,  1,  -2,  8,       0, 0 },
+                {  0, -3, -12, -1,       0, 0 },
+                {  0,  0,   0,  0, epsilon, 1 },
+                {  0,  0,   0,  0,       1, 1 }
+        }, new double[] { 2, -9, 2, 2, 1 + epsilon * epsilon, 2});
+
+        LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
+        VectorialPointValuePair optimum =
+            optimizer.optimize(problem, problem.target, new double[] { 1, 1, 1, 1, 1, 1 },
+                               new double[] { 0, 0, 0, 0, 0, 0 });
+        assertEquals(0, optimizer.getRMS(), 1.0e-10);
+        assertEquals( 3.0, optimum.getPoint()[0], 1.0e-10);
+        assertEquals( 4.0, optimum.getPoint()[1], 1.0e-10);
+        assertEquals(-1.0, optimum.getPoint()[2], 1.0e-10);
+        assertEquals(-2.0, optimum.getPoint()[3], 1.0e-10);
+        assertEquals( 1.0 + epsilon, optimum.getPoint()[4], 1.0e-10);
+        assertEquals( 1.0 - epsilon, optimum.getPoint()[5], 1.0e-10);
+
+    }
+
+    public void testNonInversible() throws ObjectiveException, OptimizationException {
+
+        LinearProblem problem = new LinearProblem(new double[][] {
+                {  1, 2, -3 },
+                {  2, 1,  3 },
+                { -3, 0, -9 }
+        }, new double[] { 1, 1, 1 });
+ 
+        LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
+        optimizer.optimize(problem, problem.target, new double[] { 1, 1, 1 }, new double[] { 0, 0, 0 });
+        assertTrue(Math.sqrt(problem.target.length) * optimizer.getRMS() > 0.6);
+        try {
+            optimizer.getCovariances();
+            fail("an exception should have been thrown");
+        } catch (OptimizationException ee) {
+            // expected behavior
+        } catch (Exception e) {
+            fail("wrong exception caught");
+        }
+
+    }
+
+    public void testIllConditioned() throws ObjectiveException, OptimizationException {
+        LinearProblem problem1 = new LinearProblem(new double[][] {
+                { 10.0, 7.0,  8.0,  7.0 },
+                {  7.0, 5.0,  6.0,  5.0 },
+                {  8.0, 6.0, 10.0,  9.0 },
+                {  7.0, 5.0,  9.0, 10.0 }
+        }, new double[] { 32, 23, 33, 31 });
+        LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
+        VectorialPointValuePair optimum1 =
+            optimizer.optimize(problem1, problem1.target, new double[] { 1, 1, 1, 1 },
+                               new double[] { 0, 1, 2, 3 });
+        assertEquals(0, optimizer.getRMS(), 1.0e-10);
+        assertEquals(1.0, optimum1.getPoint()[0], 1.0e-10);
+        assertEquals(1.0, optimum1.getPoint()[1], 1.0e-10);
+        assertEquals(1.0, optimum1.getPoint()[2], 1.0e-10);
+        assertEquals(1.0, optimum1.getPoint()[3], 1.0e-10);
+
+        LinearProblem problem2 = new LinearProblem(new double[][] {
+                { 10.00, 7.00, 8.10, 7.20 },
+                {  7.08, 5.04, 6.00, 5.00 },
+                {  8.00, 5.98, 9.89, 9.00 },
+                {  6.99, 4.99, 9.00, 9.98 }
+        }, new double[] { 32, 23, 33, 31 });
+        VectorialPointValuePair optimum2 =
+            optimizer.optimize(problem2, problem2.target, new double[] { 1, 1, 1, 1 },
+                               new double[] { 0, 1, 2, 3 });
+        assertEquals(0, optimizer.getRMS(), 1.0e-10);
+        assertEquals(-81.0, optimum2.getPoint()[0], 1.0e-8);
+        assertEquals(137.0, optimum2.getPoint()[1], 1.0e-8);
+        assertEquals(-34.0, optimum2.getPoint()[2], 1.0e-8);
+        assertEquals( 22.0, optimum2.getPoint()[3], 1.0e-8);
+
+    }
+
+    public void testMoreEstimatedParametersSimple() throws ObjectiveException, OptimizationException {
+
+        LinearProblem problem = new LinearProblem(new double[][] {
+                { 3.0, 2.0,  0.0, 0.0 },
+                { 0.0, 1.0, -1.0, 1.0 },
+                { 2.0, 0.0,  1.0, 0.0 }
+        }, new double[] { 7.0, 3.0, 5.0 });
+
+        LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
+        optimizer.optimize(problem, problem.target, new double[] { 1, 1, 1 },
+                new double[] { 7, 6, 5, 4 });
+        assertEquals(0, optimizer.getRMS(), 1.0e-10);
+
+    }
+
+    public void testMoreEstimatedParametersUnsorted() throws ObjectiveException, OptimizationException {
+        LinearProblem problem = new LinearProblem(new double[][] {
+                { 1.0, 1.0,  0.0,  0.0, 0.0,  0.0 },
+                { 0.0, 0.0,  1.0,  1.0, 1.0,  0.0 },
+                { 0.0, 0.0,  0.0,  0.0, 1.0, -1.0 },
+                { 0.0, 0.0, -1.0,  1.0, 0.0,  1.0 },
+                { 0.0, 0.0,  0.0, -1.0, 1.0,  0.0 }
+       }, new double[] { 3.0, 12.0, -1.0, 7.0, 1.0 });
+
+        LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
+        VectorialPointValuePair optimum =
+            optimizer.optimize(problem, problem.target, new double[] { 1, 1, 1, 1, 1 },
+                               new double[] { 2, 2, 2, 2, 2, 2 });
+        assertEquals(0, optimizer.getRMS(), 1.0e-10);
+        assertEquals(3.0, optimum.getPointRef()[2], 1.0e-10);
+        assertEquals(4.0, optimum.getPointRef()[3], 1.0e-10);
+        assertEquals(5.0, optimum.getPointRef()[4], 1.0e-10);
+        assertEquals(6.0, optimum.getPointRef()[5], 1.0e-10);
+
+    }
+
+    public void testRedundantEquations() throws ObjectiveException, OptimizationException {
+        LinearProblem problem = new LinearProblem(new double[][] {
+                { 1.0,  1.0 },
+                { 1.0, -1.0 },
+                { 1.0,  3.0 }
+        }, new double[] { 3.0, 1.0, 5.0 });
+
+        LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
+        VectorialPointValuePair optimum =
+            optimizer.optimize(problem, problem.target, new double[] { 1, 1, 1 },
+                               new double[] { 1, 1 });
+        assertEquals(0, optimizer.getRMS(), 1.0e-10);
+        assertEquals(2.0, optimum.getPointRef()[0], 1.0e-10);
+        assertEquals(1.0, optimum.getPointRef()[1], 1.0e-10);
+
+    }
+
+    public void testInconsistentEquations() throws ObjectiveException, OptimizationException {
+        LinearProblem problem = new LinearProblem(new double[][] {
+                { 1.0,  1.0 },
+                { 1.0, -1.0 },
+                { 1.0,  3.0 }
+        }, new double[] { 3.0, 1.0, 4.0 });
+
+        LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
+        optimizer.optimize(problem, problem.target, new double[] { 1, 1, 1 }, new double[] { 1, 1 });
+        assertTrue(optimizer.getRMS() > 0.1);
+
+    }
+
+    public void testInconsistentSizes() throws ObjectiveException, OptimizationException {
+        LinearProblem problem =
+            new LinearProblem(new double[][] { { 1, 0 }, { 0, 1 } }, new double[] { -1, 1 });
+        LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
+
+        VectorialPointValuePair optimum =
+            optimizer.optimize(problem, problem.target, new double[] { 1, 1 }, new double[] { 0, 0 });
+        assertEquals(0, optimizer.getRMS(), 1.0e-10);
+        assertEquals(-1, optimum.getPoint()[0], 1.0e-10);
+        assertEquals(+1, optimum.getPoint()[1], 1.0e-10);
+
+        try {
+            optimizer.optimize(problem, problem.target,
+                               new double[] { 1 },
+                               new double[] { 0, 0 });
+            fail("an exception should have been thrown");
+        } catch (OptimizationException oe) {
+            // expected behavior
+        } catch (Exception e) {
+            fail("wrong exception caught");
+        }
+
+        try {
+            optimizer.optimize(problem, new double[] { 1 },
+                               new double[] { 1 },
+                               new double[] { 0, 0 });
+            fail("an exception should have been thrown");
+        } catch (ObjectiveException oe) {
+            // expected behavior
+        } catch (Exception e) {
+            fail("wrong exception caught");
+        }
+
+    }
+
+    public void testControlParameters() throws OptimizationException {
+        Circle circle = new Circle();
+        circle.addPoint( 30.0,  68.0);
+        circle.addPoint( 50.0,  -6.0);
+        circle.addPoint(110.0, -20.0);
+        circle.addPoint( 35.0,  15.0);
+        circle.addPoint( 45.0,  97.0);
+        checkEstimate(circle, 0.1, 10, 1.0e-14, 1.0e-16, 1.0e-10, false);
+        checkEstimate(circle, 0.1, 10, 1.0e-15, 1.0e-17, 1.0e-10, true);
+        checkEstimate(circle, 0.1,  5, 1.0e-15, 1.0e-16, 1.0e-10, true);
+        circle.addPoint(300, -300);
+        checkEstimate(circle, 0.1, 20, 1.0e-18, 1.0e-16, 1.0e-10, true);
+    }
+
+    private void checkEstimate(VectorialDifferentiableObjectiveFunction problem,
+                               double initialStepBoundFactor, int maxCostEval,
+                               double costRelativeTolerance, double parRelativeTolerance,
+                               double orthoTolerance, boolean shouldFail) {
+        try {
+            LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
+            optimizer.setInitialStepBoundFactor(initialStepBoundFactor);
+            optimizer.setMaxEvaluations(maxCostEval);
+            optimizer.setCostRelativeTolerance(costRelativeTolerance);
+            optimizer.setParRelativeTolerance(parRelativeTolerance);
+            optimizer.setOrthoTolerance(orthoTolerance);
+            optimizer.optimize(problem, new double[] { 0, 0, 0, 0, 0 }, new double[] { 1, 1, 1, 1, 1 },
+                               new double[] { 98.680, 47.345 });
+            assertTrue(! shouldFail);
+        } catch (OptimizationException ee) {
+            assertTrue(shouldFail);
+        } catch (ObjectiveException ee) {
+            assertTrue(shouldFail);
+        } catch (Exception e) {
+            fail("wrong exception type caught");
+        }
+    }
+
+    public void testCircleFitting() throws ObjectiveException, OptimizationException {
+        Circle circle = new Circle();
+        circle.addPoint( 30.0,  68.0);
+        circle.addPoint( 50.0,  -6.0);
+        circle.addPoint(110.0, -20.0);
+        circle.addPoint( 35.0,  15.0);
+        circle.addPoint( 45.0,  97.0);
+        LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
+        VectorialPointValuePair optimum =
+            optimizer.optimize(circle, new double[] { 0, 0, 0, 0, 0 }, new double[] { 1, 1, 1, 1, 1 },
+                               new double[] { 98.680, 47.345 });
+        assertTrue(optimizer.getEvaluations() < 10);
+        assertTrue(optimizer.getJacobianEvaluations() < 10);
+        double rms = optimizer.getRMS();
+        assertEquals(1.768262623567235,  Math.sqrt(circle.getN()) * rms,  1.0e-10);
+        Point2D.Double center = new Point2D.Double(optimum.getPointRef()[0], optimum.getPointRef()[1]);
+        assertEquals(69.96016176931406, circle.getRadius(center), 1.0e-10);
+        assertEquals(96.07590211815305, center.x,      1.0e-10);
+        assertEquals(48.13516790438953, center.y,      1.0e-10);
+        double[][] cov = optimizer.getCovariances();
+        assertEquals(1.839, cov[0][0], 0.001);
+        assertEquals(0.731, cov[0][1], 0.001);
+        assertEquals(cov[0][1], cov[1][0], 1.0e-14);
+        assertEquals(0.786, cov[1][1], 0.001);
+        double[] errors = optimizer.guessParametersErrors();
+        assertEquals(1.384, errors[0], 0.001);
+        assertEquals(0.905, errors[1], 0.001);
+
+        // add perfect measurements and check errors are reduced
+        double  r = circle.getRadius(center);
+        for (double d= 0; d < 2 * Math.PI; d += 0.01) {
+            circle.addPoint(center.x + r * Math.cos(d), center.y + r * Math.sin(d));
+        }
+        double[] target = new double[circle.getN()];
+        Arrays.fill(target, 0.0);
+        double[] weights = new double[circle.getN()];
+        Arrays.fill(weights, 2.0);
+        optimum =
+            optimizer.optimize(circle, target, weights, new double[] { 98.680, 47.345 });
+        cov = optimizer.getCovariances();
+        assertEquals(0.0016, cov[0][0], 0.001);
+        assertEquals(3.2e-7, cov[0][1], 1.0e-9);
+        assertEquals(cov[0][1], cov[1][0], 1.0e-14);
+        assertEquals(0.0016, cov[1][1], 0.001);
+        errors = optimizer.guessParametersErrors();
+        assertEquals(0.002, errors[0], 0.001);
+        assertEquals(0.002, errors[1], 0.001);
+
+    }
+
+    public void testCircleFittingBadInit() throws ObjectiveException, OptimizationException {
+        Circle circle = new Circle();
+        double[][] points = new double[][] {
+                {-0.312967,  0.072366}, {-0.339248,  0.132965}, {-0.379780,  0.202724},
+                {-0.390426,  0.260487}, {-0.361212,  0.328325}, {-0.346039,  0.392619},
+                {-0.280579,  0.444306}, {-0.216035,  0.470009}, {-0.149127,  0.493832},
+                {-0.075133,  0.483271}, {-0.007759,  0.452680}, { 0.060071,  0.410235},
+                { 0.103037,  0.341076}, { 0.118438,  0.273884}, { 0.131293,  0.192201},
+                { 0.115869,  0.129797}, { 0.072223,  0.058396}, { 0.022884,  0.000718},
+                {-0.053355, -0.020405}, {-0.123584, -0.032451}, {-0.216248, -0.032862},
+                {-0.278592, -0.005008}, {-0.337655,  0.056658}, {-0.385899,  0.112526},
+                {-0.405517,  0.186957}, {-0.415374,  0.262071}, {-0.387482,  0.343398},
+                {-0.347322,  0.397943}, {-0.287623,  0.458425}, {-0.223502,  0.475513},
+                {-0.135352,  0.478186}, {-0.061221,  0.483371}, { 0.003711,  0.422737},
+                { 0.065054,  0.375830}, { 0.108108,  0.297099}, { 0.123882,  0.222850},
+                { 0.117729,  0.134382}, { 0.085195,  0.056820}, { 0.029800, -0.019138},
+                {-0.027520, -0.072374}, {-0.102268, -0.091555}, {-0.200299, -0.106578},
+                {-0.292731, -0.091473}, {-0.356288, -0.051108}, {-0.420561,  0.014926},
+                {-0.471036,  0.074716}, {-0.488638,  0.182508}, {-0.485990,  0.254068},
+                {-0.463943,  0.338438}, {-0.406453,  0.404704}, {-0.334287,  0.466119},
+                {-0.254244,  0.503188}, {-0.161548,  0.495769}, {-0.075733,  0.495560},
+                { 0.001375,  0.434937}, { 0.082787,  0.385806}, { 0.115490,  0.323807},
+                { 0.141089,  0.223450}, { 0.138693,  0.131703}, { 0.126415,  0.049174},
+                { 0.066518, -0.010217}, {-0.005184, -0.070647}, {-0.080985, -0.103635},
+                {-0.177377, -0.116887}, {-0.260628, -0.100258}, {-0.335756, -0.056251},
+                {-0.405195, -0.000895}, {-0.444937,  0.085456}, {-0.484357,  0.175597},
+                {-0.472453,  0.248681}, {-0.438580,  0.347463}, {-0.402304,  0.422428},
+                {-0.326777,  0.479438}, {-0.247797,  0.505581}, {-0.152676,  0.519380},
+                {-0.071754,  0.516264}, { 0.015942,  0.472802}, { 0.076608,  0.419077},
+                { 0.127673,  0.330264}, { 0.159951,  0.262150}, { 0.153530,  0.172681},
+                { 0.140653,  0.089229}, { 0.078666,  0.024981}, { 0.023807, -0.037022},
+                {-0.048837, -0.077056}, {-0.127729, -0.075338}, {-0.221271, -0.067526}
+        };
+        double[] target = new double[points.length];
+        Arrays.fill(target, 0.0);
+        double[] weights = new double[points.length];
+        Arrays.fill(weights, 2.0);
+        for (int i = 0; i < points.length; ++i) {
+            circle.addPoint(points[i][0], points[i][1]);
+        }
+        LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
+        optimizer.setConvergenceChecker(new SimpleVectorialValueChecker(1.0e-10, 1.0e-10));
+        VectorialPointValuePair optimum =
+            optimizer.optimize(circle, target, weights, new double[] { -12, -12 });
+        Point2D.Double center = new Point2D.Double(optimum.getPointRef()[0], optimum.getPointRef()[1]);
+        assertTrue(optimizer.getEvaluations() < 25);
+        assertTrue(optimizer.getJacobianEvaluations() < 20);
+        assertEquals( 0.043, optimizer.getRMS(), 1.0e-3);
+        assertEquals( 0.292235,  circle.getRadius(center), 1.0e-6);
+        assertEquals(-0.151738,  center.x,      1.0e-6);
+        assertEquals( 0.2075001, center.y,      1.0e-6);
+    }
+
+    public void testMath199() throws ObjectiveException, OptimizationException {
+        try {
+            QuadraticProblem problem = new QuadraticProblem();
+            problem.addPoint (0, -3.182591015485607);
+            problem.addPoint (1, -2.5581184967730577);
+            problem.addPoint (2, -2.1488478161387325);
+            problem.addPoint (3, -1.9122489313410047);
+            problem.addPoint (4, 1.7785661310051026);
+            new LevenbergMarquardtOptimizer().optimize(problem,
+                                                       new double[] { 0, 0, 0, 0, 0 },
+                                                       new double[] { 0.0, 4.4e-323, 1.0, 4.4e-323, 0.0 },
+                                                       new double[] { 0, 0, 0 });
+            fail("an exception should have been thrown");
+        } catch (OptimizationException ee) {
+            // expected behavior
+        }
+
+    }
 
-  public void testOneSet() throws OptimizationException {
+    private static class LinearProblem implements VectorialDifferentiableObjectiveFunction {
 
-    EstimatedParameter[] p = {
-       new EstimatedParameter("p0", 0),
-       new EstimatedParameter("p1", 0),
-       new EstimatedParameter("p2", 0)
-    };
-    LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
-      new LinearMeasurement(new double[] { 1.0 },
-                            new EstimatedParameter[] { p[0] },
-                            1.0),
-      new LinearMeasurement(new double[] { -1.0, 1.0 },
-                            new EstimatedParameter[] { p[0], p[1] },
-                            1.0),
-      new LinearMeasurement(new double[] { -1.0, 1.0 },
-                            new EstimatedParameter[] { p[1], p[2] },
-                            1.0)
-    });
-
-    LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
-    estimator.estimate(problem);
-    assertEquals(0, estimator.getRMS(problem), 1.0e-10);
-    assertEquals(1.0, p[0].getEstimate(), 1.0e-10);
-    assertEquals(2.0, p[1].getEstimate(), 1.0e-10);
-    assertEquals(3.0, p[2].getEstimate(), 1.0e-10);
-
-  }
-
-  public void testTwoSets() throws OptimizationException {
-    EstimatedParameter[] p = {
-      new EstimatedParameter("p0", 0),
-      new EstimatedParameter("p1", 1),
-      new EstimatedParameter("p2", 2),
-      new EstimatedParameter("p3", 3),
-      new EstimatedParameter("p4", 4),
-      new EstimatedParameter("p5", 5)
-    };
-
-    double epsilon = 1.0e-7;
-    LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
-
-      // 4 elements sub-problem
-      new LinearMeasurement(new double[] {  2.0,  1.0,  4.0 },
-                            new EstimatedParameter[] { p[0], p[1], p[3] },
-                            2.0),
-      new LinearMeasurement(new double[] { -4.0, -2.0,   3.0, -7.0 },
-                           new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
-                           -9.0),
-      new LinearMeasurement(new double[] {  4.0,  1.0,  -2.0,  8.0 },
-                            new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
-                            2.0),
-      new LinearMeasurement(new double[] { -3.0, -12.0, -1.0 },
-                           new EstimatedParameter[] { p[1], p[2], p[3] },
-                           2.0),
-
-      // 2 elements sub-problem
-      new LinearMeasurement(new double[] { epsilon, 1.0 },
-                            new EstimatedParameter[] { p[4], p[5] },
-                            1.0 + epsilon * epsilon),
-      new LinearMeasurement(new double[] {  1.0, 1.0 },
-                            new EstimatedParameter[] { p[4], p[5] },
-                            2.0)
-
-    });
-
-    LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
-    estimator.estimate(problem);
-    assertEquals(0, estimator.getRMS(problem), 1.0e-10);
-    assertEquals( 3.0, p[0].getEstimate(), 1.0e-10);
-    assertEquals( 4.0, p[1].getEstimate(), 1.0e-10);
-    assertEquals(-1.0, p[2].getEstimate(), 1.0e-10);
-    assertEquals(-2.0, p[3].getEstimate(), 1.0e-10);
-    assertEquals( 1.0 + epsilon, p[4].getEstimate(), 1.0e-10);
-    assertEquals( 1.0 - epsilon, p[5].getEstimate(), 1.0e-10);
-
-  }
-
-  public void testNonInversible() throws OptimizationException {
-
-    EstimatedParameter[] p = {
-       new EstimatedParameter("p0", 0),
-       new EstimatedParameter("p1", 0),
-       new EstimatedParameter("p2", 0)
-    };
-    LinearMeasurement[] m = new LinearMeasurement[] {
-      new LinearMeasurement(new double[] {  1.0, 2.0, -3.0 },
-                            new EstimatedParameter[] { p[0], p[1], p[2] },
-                            1.0),
-      new LinearMeasurement(new double[] {  2.0, 1.0,  3.0 },
-                            new EstimatedParameter[] { p[0], p[1], p[2] },
-                            1.0),
-      new LinearMeasurement(new double[] { -3.0, -9.0 },
-                            new EstimatedParameter[] { p[0], p[2] },
-                            1.0)
-    };
-    LinearProblem problem = new LinearProblem(m);
-
-    LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
-    double initialCost = estimator.getRMS(problem);
-    estimator.estimate(problem);
-    assertTrue(estimator.getRMS(problem) < initialCost);
-    assertTrue(Math.sqrt(m.length) * estimator.getRMS(problem) > 0.6);
-    try {
-        estimator.getCovariances(problem);
-        fail("an exception should have been thrown");
-    } catch (OptimizationException ee) {
-        // expected behavior
-    } catch (Exception e) {
-        fail("wrong exception caught");
-    }
-   double dJ0 = 2 * (m[0].getResidual() * m[0].getPartial(p[0])
-                    + m[1].getResidual() * m[1].getPartial(p[0])
-                    + m[2].getResidual() * m[2].getPartial(p[0]));
-    double dJ1 = 2 * (m[0].getResidual() * m[0].getPartial(p[1])
-                    + m[1].getResidual() * m[1].getPartial(p[1]));
-    double dJ2 = 2 * (m[0].getResidual() * m[0].getPartial(p[2])
-                    + m[1].getResidual() * m[1].getPartial(p[2])
-                    + m[2].getResidual() * m[2].getPartial(p[2]));
-    assertEquals(0, dJ0, 1.0e-10);
-    assertEquals(0, dJ1, 1.0e-10);
-    assertEquals(0, dJ2, 1.0e-10);
-
-  }
-
-  public void testIllConditioned() throws OptimizationException {
-    EstimatedParameter[] p = {
-      new EstimatedParameter("p0", 0),
-      new EstimatedParameter("p1", 1),
-      new EstimatedParameter("p2", 2),
-      new EstimatedParameter("p3", 3)
-    };
-
-    LinearProblem problem1 = new LinearProblem(new LinearMeasurement[] {
-      new LinearMeasurement(new double[] { 10.0, 7.0,  8.0,  7.0 },
-                            new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
-                            32.0),
-      new LinearMeasurement(new double[] {  7.0, 5.0,  6.0,  5.0 },
-                            new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
-                            23.0),
-      new LinearMeasurement(new double[] {  8.0, 6.0, 10.0,  9.0 },
-                            new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
-                            33.0),
-      new LinearMeasurement(new double[] {  7.0, 5.0,  9.0, 10.0 },
-                            new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
-                            31.0)
-    });
-    LevenbergMarquardtEstimator estimator1 = new LevenbergMarquardtEstimator();
-    estimator1.estimate(problem1);
-    assertEquals(0, estimator1.getRMS(problem1), 1.0e-10);
-    assertEquals(1.0, p[0].getEstimate(), 1.0e-10);
-    assertEquals(1.0, p[1].getEstimate(), 1.0e-10);
-    assertEquals(1.0, p[2].getEstimate(), 1.0e-10);
-    assertEquals(1.0, p[3].getEstimate(), 1.0e-10);
-
-    LinearProblem problem2 = new LinearProblem(new LinearMeasurement[] {
-      new LinearMeasurement(new double[] { 10.0, 7.0,  8.1,  7.2 },
-                            new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
-                            32.0),
-      new LinearMeasurement(new double[] {  7.08, 5.04,  6.0,  5.0 },
-                            new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
-                            23.0),
-      new LinearMeasurement(new double[] {  8.0, 5.98, 9.89,  9.0 },
-                             new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
-                            33.0),
-      new LinearMeasurement(new double[] {  6.99, 4.99,  9.0, 9.98 },
-                             new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
-                            31.0)
-    });
-    LevenbergMarquardtEstimator estimator2 = new LevenbergMarquardtEstimator();
-    estimator2.estimate(problem2);
-    assertEquals(0, estimator2.getRMS(problem2), 1.0e-10);
-    assertEquals(-81.0, p[0].getEstimate(), 1.0e-8);
-    assertEquals(137.0, p[1].getEstimate(), 1.0e-8);
-    assertEquals(-34.0, p[2].getEstimate(), 1.0e-8);
-    assertEquals( 22.0, p[3].getEstimate(), 1.0e-8);
-
-  }
-
-  public void testMoreEstimatedParametersSimple() throws OptimizationException {
-
-    EstimatedParameter[] p = {
-       new EstimatedParameter("p0", 7),
-       new EstimatedParameter("p1", 6),
-       new EstimatedParameter("p2", 5),
-       new EstimatedParameter("p3", 4)
-     };
-    LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
-      new LinearMeasurement(new double[] { 3.0, 2.0 },
-                             new EstimatedParameter[] { p[0], p[1] },
-                             7.0),
-      new LinearMeasurement(new double[] { 1.0, -1.0, 1.0 },
-                             new EstimatedParameter[] { p[1], p[2], p[3] },
-                             3.0),
-      new LinearMeasurement(new double[] { 2.0, 1.0 },
-                             new EstimatedParameter[] { p[0], p[2] },
-                             5.0)
-    });
-
-    LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
-    estimator.estimate(problem);
-    assertEquals(0, estimator.getRMS(problem), 1.0e-10);
-
-  }
-
-  public void testMoreEstimatedParametersUnsorted() throws OptimizationException {
-    EstimatedParameter[] p = {
-      new EstimatedParameter("p0", 2),
-      new EstimatedParameter("p1", 2),
-      new EstimatedParameter("p2", 2),
-      new EstimatedParameter("p3", 2),
-      new EstimatedParameter("p4", 2),
-      new EstimatedParameter("p5", 2)
-    };
-    LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
-      new LinearMeasurement(new double[] { 1.0, 1.0 },
-                           new EstimatedParameter[] { p[0], p[1] },
-                           3.0),
-      new LinearMeasurement(new double[] { 1.0, 1.0, 1.0 },
-                           new EstimatedParameter[] { p[2], p[3], p[4] },
-                           12.0),
-      new LinearMeasurement(new double[] { 1.0, -1.0 },
-                           new EstimatedParameter[] { p[4], p[5] },
-                           -1.0),
-      new LinearMeasurement(new double[] { 1.0, -1.0, 1.0 },
-                           new EstimatedParameter[] { p[3], p[2], p[5] },
-                           7.0),
-      new LinearMeasurement(new double[] { 1.0, -1.0 },
-                           new EstimatedParameter[] { p[4], p[3] },
-                           1.0)
-    });
-
-    LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
-    estimator.estimate(problem);
-    assertEquals(0, estimator.getRMS(problem), 1.0e-10);
-    assertEquals(3.0, p[2].getEstimate(), 1.0e-10);
-    assertEquals(4.0, p[3].getEstimate(), 1.0e-10);
-    assertEquals(5.0, p[4].getEstimate(), 1.0e-10);
-    assertEquals(6.0, p[5].getEstimate(), 1.0e-10);
-
-  }
-
-  public void testRedundantEquations() throws OptimizationException {
-    EstimatedParameter[] p = {
-      new EstimatedParameter("p0", 1),
-      new EstimatedParameter("p1", 1)
-    };
-    LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
-      new LinearMeasurement(new double[] { 1.0, 1.0 },
-                             new EstimatedParameter[] { p[0], p[1] },
-                             3.0),
-      new LinearMeasurement(new double[] { 1.0, -1.0 },
-                             new EstimatedParameter[] { p[0], p[1] },
-                             1.0),
-      new LinearMeasurement(new double[] { 1.0, 3.0 },
-                             new EstimatedParameter[] { p[0], p[1] },
-                             5.0)
-    });
-
-    LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
-    estimator.estimate(problem);
-    assertEquals(0, estimator.getRMS(problem), 1.0e-10);
-    assertEquals(2.0, p[0].getEstimate(), 1.0e-10);
-    assertEquals(1.0, p[1].getEstimate(), 1.0e-10);
-
-  }
-
-  public void testInconsistentEquations() throws OptimizationException {
-    EstimatedParameter[] p = {
-      new EstimatedParameter("p0", 1),
-      new EstimatedParameter("p1", 1)
-    };
-    LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
-      new LinearMeasurement(new double[] { 1.0, 1.0 },
-                            new EstimatedParameter[] { p[0], p[1] },
-                            3.0),
-      new LinearMeasurement(new double[] { 1.0, -1.0 },
-                            new EstimatedParameter[] { p[0], p[1] },
-                            1.0),
-      new LinearMeasurement(new double[] { 1.0, 3.0 },
-                            new EstimatedParameter[] { p[0], p[1] },
-                            4.0)
-    });
-
-    LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
-    estimator.estimate(problem);
-    assertTrue(estimator.getRMS(problem) > 0.1);
-
-  }
-
-  public void testControlParameters() throws OptimizationException {
-      Circle circle = new Circle(98.680, 47.345);
-      circle.addPoint( 30.0,  68.0);
-      circle.addPoint( 50.0,  -6.0);
-      circle.addPoint(110.0, -20.0);
-      circle.addPoint( 35.0,  15.0);
-      circle.addPoint( 45.0,  97.0);
-      checkEstimate(circle, 0.1, 10, 1.0e-14, 1.0e-16, 1.0e-10, false);
-      checkEstimate(circle, 0.1, 10, 1.0e-15, 1.0e-17, 1.0e-10, true);
-      checkEstimate(circle, 0.1,  5, 1.0e-15, 1.0e-16, 1.0e-10, true);
-      circle.addPoint(300, -300);
-      checkEstimate(circle, 0.1, 20, 1.0e-18, 1.0e-16, 1.0e-10, true);
-  }
-
-  private void checkEstimate(EstimationProblem problem,
-                             double initialStepBoundFactor, int maxCostEval,
-                             double costRelativeTolerance, double parRelativeTolerance,
-                             double orthoTolerance, boolean shouldFail) {
-      try {
-        LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
-        estimator.setInitialStepBoundFactor(initialStepBoundFactor);
-        estimator.setMaxCostEval(maxCostEval);
-        estimator.setCostRelativeTolerance(costRelativeTolerance);
-        estimator.setParRelativeTolerance(parRelativeTolerance);
-        estimator.setOrthoTolerance(orthoTolerance);
-        estimator.estimate(problem);
-        assertTrue(! shouldFail);
-      } catch (OptimizationException ee) {
-        assertTrue(shouldFail);
-      } catch (Exception e) {
-        fail("wrong exception type caught");
-      }
-    }
-
-  public void testCircleFitting() throws OptimizationException {
-      Circle circle = new Circle(98.680, 47.345);
-      circle.addPoint( 30.0,  68.0);
-      circle.addPoint( 50.0,  -6.0);
-      circle.addPoint(110.0, -20.0);
-      circle.addPoint( 35.0,  15.0);
-      circle.addPoint( 45.0,  97.0);
-      LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
-      estimator.estimate(circle);
-      assertTrue(estimator.getCostEvaluations() < 10);
-      assertTrue(estimator.getJacobianEvaluations() < 10);
-      double rms = estimator.getRMS(circle);
-      assertEquals(1.768262623567235,  Math.sqrt(circle.getM()) * rms,  1.0e-10);
-      assertEquals(69.96016176931406, circle.getRadius(), 1.0e-10);
-      assertEquals(96.07590211815305, circle.getX(),      1.0e-10);
-      assertEquals(48.13516790438953, circle.getY(),      1.0e-10);
-      double[][] cov = estimator.getCovariances(circle);
-      assertEquals(1.839, cov[0][0], 0.001);
-      assertEquals(0.731, cov[0][1], 0.001);
-      assertEquals(cov[0][1], cov[1][0], 1.0e-14);
-      assertEquals(0.786, cov[1][1], 0.001);
-      double[] errors = estimator.guessParametersErrors(circle);
-      assertEquals(1.384, errors[0], 0.001);
-      assertEquals(0.905, errors[1], 0.001);
-  
-      // add perfect measurements and check errors are reduced
-      double cx = circle.getX();
-      double cy = circle.getY();
-      double  r = circle.getRadius();
-      for (double d= 0; d < 2 * Math.PI; d += 0.01) {
-          circle.addPoint(cx + r * Math.cos(d), cy + r * Math.sin(d));
-      }
-      estimator = new LevenbergMarquardtEstimator();
-      estimator.estimate(circle);
-      cov = estimator.getCovariances(circle);
-      assertEquals(0.004, cov[0][0], 0.001);
-      assertEquals(6.40e-7, cov[0][1], 1.0e-9);
-      assertEquals(cov[0][1], cov[1][0], 1.0e-14);
-      assertEquals(0.003, cov[1][1], 0.001);
-      errors = estimator.guessParametersErrors(circle);
-      assertEquals(0.004, errors[0], 0.001);
-      assertEquals(0.004, errors[1], 0.001);
-
-  }
-
-  public void testCircleFittingBadInit() throws OptimizationException {
-    Circle circle = new Circle(-12, -12);
-    double[][] points = new double[][] {
-      {-0.312967,  0.072366}, {-0.339248,  0.132965}, {-0.379780,  0.202724},
-      {-0.390426,  0.260487}, {-0.361212,  0.328325}, {-0.346039,  0.392619},
-      {-0.280579,  0.444306}, {-0.216035,  0.470009}, {-0.149127,  0.493832},
-      {-0.075133,  0.483271}, {-0.007759,  0.452680}, { 0.060071,  0.410235},
-      { 0.103037,  0.341076}, { 0.118438,  0.273884}, { 0.131293,  0.192201},
-      { 0.115869,  0.129797}, { 0.072223,  0.058396}, { 0.022884,  0.000718},
-      {-0.053355, -0.020405}, {-0.123584, -0.032451}, {-0.216248, -0.032862},
-      {-0.278592, -0.005008}, {-0.337655,  0.056658}, {-0.385899,  0.112526},
-      {-0.405517,  0.186957}, {-0.415374,  0.262071}, {-0.387482,  0.343398},
-      {-0.347322,  0.397943}, {-0.287623,  0.458425}, {-0.223502,  0.475513},
-      {-0.135352,  0.478186}, {-0.061221,  0.483371}, { 0.003711,  0.422737},
-      { 0.065054,  0.375830}, { 0.108108,  0.297099}, { 0.123882,  0.222850},
-      { 0.117729,  0.134382}, { 0.085195,  0.056820}, { 0.029800, -0.019138},
-      {-0.027520, -0.072374}, {-0.102268, -0.091555}, {-0.200299, -0.106578},
-      {-0.292731, -0.091473}, {-0.356288, -0.051108}, {-0.420561,  0.014926},
-      {-0.471036,  0.074716}, {-0.488638,  0.182508}, {-0.485990,  0.254068},
-      {-0.463943,  0.338438}, {-0.406453,  0.404704}, {-0.334287,  0.466119},
-      {-0.254244,  0.503188}, {-0.161548,  0.495769}, {-0.075733,  0.495560},
-      { 0.001375,  0.434937}, { 0.082787,  0.385806}, { 0.115490,  0.323807},
-      { 0.141089,  0.223450}, { 0.138693,  0.131703}, { 0.126415,  0.049174},
-      { 0.066518, -0.010217}, {-0.005184, -0.070647}, {-0.080985, -0.103635},
-      {-0.177377, -0.116887}, {-0.260628, -0.100258}, {-0.335756, -0.056251},
-      {-0.405195, -0.000895}, {-0.444937,  0.085456}, {-0.484357,  0.175597},
-      {-0.472453,  0.248681}, {-0.438580,  0.347463}, {-0.402304,  0.422428},
-      {-0.326777,  0.479438}, {-0.247797,  0.505581}, {-0.152676,  0.519380},
-      {-0.071754,  0.516264}, { 0.015942,  0.472802}, { 0.076608,  0.419077},
-      { 0.127673,  0.330264}, { 0.159951,  0.262150}, { 0.153530,  0.172681},
-      { 0.140653,  0.089229}, { 0.078666,  0.024981}, { 0.023807, -0.037022},
-      {-0.048837, -0.077056}, {-0.127729, -0.075338}, {-0.221271, -0.067526}
-    };
-    for (int i = 0; i < points.length; ++i) {
-      circle.addPoint(points[i][0], points[i][1]);
-    }
-    LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
-    estimator.estimate(circle);
-    assertTrue(estimator.getCostEvaluations() < 15);
-    assertTrue(estimator.getJacobianEvaluations() < 10);
-    assertEquals( 0.030184491196225207, estimator.getRMS(circle), 1.0e-9);
-    assertEquals( 0.2922350065939634,   circle.getRadius(), 1.0e-9);
-    assertEquals(-0.15173845023862165,  circle.getX(),      1.0e-8);
-    assertEquals( 0.20750021499570379,  circle.getY(),      1.0e-8);
-  }
-
-  public void testMath199() {
-      try {
-          QuadraticProblem problem = new QuadraticProblem();
-          problem.addPoint (0, -3.182591015485607, 0.0);
-          problem.addPoint (1, -2.5581184967730577, 4.4E-323);
-          problem.addPoint (2, -2.1488478161387325, 1.0);
-          problem.addPoint (3, -1.9122489313410047, 4.4E-323);
-          problem.addPoint (4, 1.7785661310051026, 0.0);
-          new LevenbergMarquardtEstimator().estimate(problem);
-          fail("an exception should have been thrown");
-      } catch (OptimizationException ee) {
-          // expected behavior
-      }
-
-  }
-
-  private static class LinearProblem implements EstimationProblem {
-
-    public LinearProblem(LinearMeasurement[] measurements) {
-      this.measurements = measurements;
-    }
-
-    public WeightedMeasurement[] getMeasurements() {
-      return measurements;
-    }
-
-    public EstimatedParameter[] getUnboundParameters() {
-      return getAllParameters();
-    }
-
-    public EstimatedParameter[] getAllParameters() {
-      HashSet<EstimatedParameter> set = new HashSet<EstimatedParameter>();
-      for (int i = 0; i < measurements.length; ++i) {
-        EstimatedParameter[] parameters = measurements[i].getParameters();
-        for (int j = 0; j < parameters.length; ++j) {
-          set.add(parameters[j]);
-        }
-      }
-      return (EstimatedParameter[]) set.toArray(new EstimatedParameter[set.size()]);
-    }
-  
-    private LinearMeasurement[] measurements;
-
-  }
-
-  private static class LinearMeasurement extends WeightedMeasurement {
-
-    public LinearMeasurement(double[] factors, EstimatedParameter[] parameters,
-                             double setPoint) {
-      super(1.0, setPoint);
-      this.factors = factors;
-      this.parameters = parameters;
-    }
-
-    public double getTheoreticalValue() {
-      double v = 0;
-      for (int i = 0; i < factors.length; ++i) {
-        v += factors[i] * parameters[i].getEstimate();
-      }
-      return v;
-    }
-
-    public double getPartial(EstimatedParameter parameter) {
-      for (int i = 0; i < parameters.length; ++i) {
-        if (parameters[i] == parameter) {
-          return factors[i];
-        }
-      }
-      return 0;
-    }
-
-    public EstimatedParameter[] getParameters() {
-      return parameters;
-    }
-
-    private double[] factors;
-    private EstimatedParameter[] parameters;
-    private static final long serialVersionUID = -3922448707008868580L;
-
-  }
-
-  private static class Circle implements EstimationProblem {
-
-    public Circle(double cx, double cy) {
-      this.cx = new EstimatedParameter("cx", cx);
-      this.cy = new EstimatedParameter("cy", cy);
-      points  = new ArrayList<PointModel>();
-    }
-
-    public void addPoint(double px, double py) {
-      points.add(new PointModel(px, py));
-    }
-
-    public int getM() {
-      return points.size();
-    }
-
-    public WeightedMeasurement[] getMeasurements() {
-      return (WeightedMeasurement[]) points.toArray(new PointModel[points.size()]);
-    }
-
-    public EstimatedParameter[] getAllParameters() {
-      return new EstimatedParameter[] { cx, cy };
-    }
-
-    public EstimatedParameter[] getUnboundParameters() {
-      return new EstimatedParameter[] { cx, cy };
-    }
-
-    public double getPartialRadiusX() {
-      double dRdX = 0;
-      for (PointModel point : points) {
-        dRdX += point.getPartialDiX();
-      }
-      return dRdX / points.size();
-    }
-
-    public double getPartialRadiusY() {
-      double dRdY = 0;
-      for (PointModel point : points) {
-        dRdY += point.getPartialDiY();
-      }
-      return dRdY / points.size();
-    }
-
-   public double getRadius() {
-      double r = 0;
-      for (PointModel point : points) {
-        r += point.getCenterDistance();
-      }
-      return r / points.size();
-    }
-
-    public double getX() {
-      return cx.getEstimate();
-    }
-
-    public double getY() {
-      return cy.getEstimate();
-    }
-
-    private class PointModel extends WeightedMeasurement {
-
-      public PointModel(double px, double py) {
-        super(1.0, 0.0);
-        this.px = px;
-        this.py = py;
-      }
-
-      public double getPartial(EstimatedParameter parameter) {
-        if (parameter == cx) {
-          return getPartialDiX() - getPartialRadiusX();
-        } else if (parameter == cy) {
-          return getPartialDiY() - getPartialRadiusY();
-        }
-        return 0;
-      }
-
-      public double getCenterDistance() {
-        double dx = px - cx.getEstimate();
-        double dy = py - cy.getEstimate();
-        return Math.sqrt(dx * dx + dy * dy);
-      }
-
-      public double getPartialDiX() {
-        return (cx.getEstimate() - px) / getCenterDistance();
-      }
-
-      public double getPartialDiY() {
-        return (cy.getEstimate() - py) / getCenterDistance();
-      }
-
-      public double getTheoreticalValue() {
-        return getCenterDistance() - getRadius();
-      }
-
-      private double px;
-      private double py;
-      private static final long serialVersionUID = 1L;
-
-    }
-
-    private EstimatedParameter cx;
-    private EstimatedParameter cy;
-    private ArrayList<PointModel> points;
-
-  }
-
-  private static class QuadraticProblem extends SimpleEstimationProblem {
-
-      private EstimatedParameter a;
-      private EstimatedParameter b;
-      private EstimatedParameter c;
-
-      public QuadraticProblem() {
-          a = new EstimatedParameter("a", 0.0);
-          b = new EstimatedParameter("b", 0.0);
-          c = new EstimatedParameter("c", 0.0);
-          addParameter(a);
-          addParameter(b);
-          addParameter(c);
-      }
-
-      public void addPoint(double x, double y, double w) {
-          addMeasurement(new LocalMeasurement(x, y, w));
-      }
-
-      public double getA() {
-          return a.getEstimate();
-      }
-
-      public double getB() {
-          return b.getEstimate();
-      }
-
-      public double getC() {
-          return c.getEstimate();
-      }
-
-      public double theoreticalValue(double x) {
-          return ( (a.getEstimate() * x + b.getEstimate() ) * x + c.getEstimate());
-      }
-
-      private double partial(double x, EstimatedParameter parameter) {
-          if (parameter == a) {
-              return x * x;
-          } else if (parameter == b) {
-              return x;
-          } else {
-              return 1.0;
-          }
-      }
-
-      private class LocalMeasurement extends WeightedMeasurement {
-
-        private static final long serialVersionUID = 1555043155023729130L;
-        private final double x;
-
-          // constructor
-          public LocalMeasurement(double x, double y, double w) {
-              super(w, y);
-              this.x = x;
-          }
-
-          public double getTheoreticalValue() {
-              return theoreticalValue(x);
-          }
-
-          public double getPartial(EstimatedParameter parameter) {
-              return partial(x, parameter);
-          }
-
-      }
-  }
-
-  public static Test suite() {
-    return new TestSuite(LevenbergMarquardtEstimatorTest.class);
-  }
+        private static final long serialVersionUID = 703247177355019415L;
+        final RealMatrix factors;
+        final double[] target;
+        public LinearProblem(double[][] factors, double[] target) {
+            this.factors = new DenseRealMatrix(factors);
+            this.target  = target;
+        }
+
+        public double[][] jacobian(double[] variables, double[] value) {
+            return factors.getData();
+        }
+
+        public double[] objective(double[] variables) {
+            return factors.operate(variables);
+        }
+
+    }
+
+    private static class Circle implements VectorialDifferentiableObjectiveFunction {
+
+        private static final long serialVersionUID = -4711170319243817874L;
+
+        private ArrayList<Point2D.Double> points;
+
+        public Circle() {
+            points  = new ArrayList<Point2D.Double>();
+        }
+
+        public void addPoint(double px, double py) {
+            points.add(new Point2D.Double(px, py));
+        }
+
+        public int getN() {
+            return points.size();
+        }
+
+        public double getRadius(Point2D.Double center) {
+            double r = 0;
+            for (Point2D.Double point : points) {
+                r += point.distance(center);
+            }
+            return r / points.size();
+        }
+
+        public double[][] jacobian(double[] variables, double[] value)
+        throws ObjectiveException, IllegalArgumentException {
+
+            int n = points.size();
+            Point2D.Double center = new Point2D.Double(variables[0], variables[1]);
+
+            // gradient of the optimal radius
+            double dRdX = 0;
+            double dRdY = 0;
+            for (Point2D.Double pk : points) {
+                double dk = pk.distance(center);
+                dRdX += (center.x - pk.x) / dk;
+                dRdY += (center.y - pk.y) / dk;
+            }
+            dRdX /= n;
+            dRdY /= n;
+
+            // jacobian of the radius residuals
+            double[][] jacobian = new double[n][2];
+            for (int i = 0; i < n; ++i) {
+                Point2D.Double pi = points.get(i);
+                double di   = pi.distance(center);
+                jacobian[i][0] = (center.x - pi.x) / di - dRdX;    
+                jacobian[i][1] = (center.y - pi.y) / di - dRdY;    
+            }
+
+            return jacobian;
+
+        }
+
+        public double[] objective(double[] variables)
+        throws ObjectiveException, IllegalArgumentException {
+
+            Point2D.Double center = new Point2D.Double(variables[0], variables[1]);
+            double radius = getRadius(center);
+
+            double[] residuals = new double[points.size()];
+            for (int i = 0; i < residuals.length; ++i) {
+                residuals[i] = points.get(i).distance(center) - radius;
+            }
+
+            return residuals;
+
+        }
+
+    }
+
+    private static class QuadraticProblem implements VectorialDifferentiableObjectiveFunction {
+
+        private static final long serialVersionUID = -247096133023967957L;
+        private List<Double> x;
+        private List<Double> y;
+
+        public QuadraticProblem() {
+            x = new ArrayList<Double>();
+            y = new ArrayList<Double>();
+        }
+
+        public void addPoint(double x, double y) {
+            this.x.add(x);
+            this.y.add(y);
+        }
+
+        public double[][] jacobian(double[] variables, double[] value) {
+            double[][] jacobian = new double[x.size()][3];
+            for (int i = 0; i < jacobian.length; ++i) {
+                jacobian[i][0] = x.get(i) * x.get(i);
+                jacobian[i][1] = x.get(i);
+                jacobian[i][2] = 1.0;
+            }
+            return jacobian;
+        }
+
+        public double[] objective(double[] variables) {
+            double[] values = new double[x.size()];
+            for (int i = 0; i < values.length; ++i) {
+                values[i] = (variables[0] * x.get(i) + variables[1]) * x.get(i) + variables[2];
+            }
+            return values;
+        }
+
+    }
+
+    public static Test suite() {
+        return new TestSuite(LevenbergMarquardtOptimizerTest.class);
+    }
 
 }



Mime
View raw message