commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From er...@apache.org
Subject svn commit: r1420684 [13/15] - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math3/exception/ main/java/org/apache/commons/math3/exception/util/ main/java/org/apache/commons/math3/fitting/ main/java/org/apache/commons/math3/optim/ mai...
Date Wed, 12 Dec 2012 14:11:04 GMT
Added: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optim/nonlinear/vector/jacobian/AbstractLeastSquaresOptimizerTestValidation.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optim/nonlinear/vector/jacobian/AbstractLeastSquaresOptimizerTestValidation.java?rev=1420684&view=auto
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optim/nonlinear/vector/jacobian/AbstractLeastSquaresOptimizerTestValidation.java (added)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optim/nonlinear/vector/jacobian/AbstractLeastSquaresOptimizerTestValidation.java Wed Dec 12 14:10:38 2012
@@ -0,0 +1,334 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with this
+ * work for additional information regarding copyright ownership. The ASF
+ * licenses this file to You under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ * http://www.apache.org/licenses/LICENSE-2.0 Unless required by applicable law
+ * or agreed to in writing, software distributed under the License is
+ * distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied. See the License for the specific language
+ * governing permissions and limitations under the License.
+ */
+package org.apache.commons.math3.optim.nonlinear.vector.jacobian;
+
+import java.util.Arrays;
+import java.util.List;
+import java.util.ArrayList;
+import java.awt.geom.Point2D;
+import org.apache.commons.math3.optim.PointVectorValuePair;
+import org.apache.commons.math3.optim.InitialGuess;
+import org.apache.commons.math3.optim.MaxEval;
+import org.apache.commons.math3.optim.nonlinear.vector.Target;
+import org.apache.commons.math3.optim.nonlinear.vector.Weight;
+import org.apache.commons.math3.stat.descriptive.SummaryStatistics;
+import org.apache.commons.math3.stat.descriptive.StatisticalSummary;
+import org.apache.commons.math3.util.FastMath;
+import org.junit.Test;
+import org.junit.Assert;
+
+/**
+ * This class demonstrates the main functionality of the
+ * {@link AbstractLeastSquaresOptimizer}, common to the
+ * optimizer implementations in package
+ * {@link org.apache.commons.math3.optimization.general}.
+ * <br/>
+ * Not enabled by default, as the class name does not end with "Test".
+ * <br/>
+ * Invoke by running
+ * <pre><code>
+ *  mvn test -Dtest=AbstractLeastSquaresOptimizerTestValidation
+ * </code></pre>
+ * or by running
+ * <pre><code>
+ *  mvn test -Dtest=AbstractLeastSquaresOptimizerTestValidation -DargLine="-DmcRuns=1234 -server"
+ * </code></pre>
+ */
+public class AbstractLeastSquaresOptimizerTestValidation {
+    private static final int MONTE_CARLO_RUNS = Integer.parseInt(System.getProperty("mcRuns",
+                                                                                    "100"));
+
+    /**
+     * Using a Monte-Carlo procedure, this test checks the error estimations
+     * as provided by the square-root of the diagonal elements of the
+     * covariance matrix.
+     * <br/>
+     * The test generates sets of observations, each sampled from
+     * a Gaussian distribution.
+     * <br/>
+     * The optimization problem solved is defined in class
+     * {@link StraightLineProblem}.
+     * <br/>
+     * The output (on stdout) will be a table summarizing the distribution
+     * of parameters generated by the Monte-Carlo process and by the direct
+     * estimation provided by the diagonal elements of the covariance matrix.
+     */
+    @Test
+    public void testParametersErrorMonteCarloObservations() {
+        // Error on the observations.
+        final double yError = 15;
+
+        // True values of the parameters.
+        final double slope = 123.456;
+        final double offset = -98.765;
+
+        // Samples generator.
+        final RandomStraightLinePointGenerator lineGenerator
+            = new RandomStraightLinePointGenerator(slope, offset,
+                                                   yError,
+                                                   -1e3, 1e4,
+                                                   138577L);
+
+        // Number of observations.
+        final int numObs = 100; // XXX Should be a command-line option.
+        // number of parameters.
+        final int numParams = 2;
+
+        // Parameters found for each of Monte-Carlo run.
+        final SummaryStatistics[] paramsFoundByDirectSolution = new SummaryStatistics[numParams];
+        // Sigma estimations (square-root of the diagonal elements of the
+        // covariance matrix), for each Monte-Carlo run.
+        final SummaryStatistics[] sigmaEstimate = new SummaryStatistics[numParams];
+
+        // Initialize statistics accumulators.
+        for (int i = 0; i < numParams; i++) {
+            paramsFoundByDirectSolution[i] = new SummaryStatistics();
+            sigmaEstimate[i] = new SummaryStatistics();
+        }
+
+        // Dummy optimizer (to compute the covariance matrix).
+        final AbstractLeastSquaresOptimizer optim = new DummyOptimizer();
+        final double[] init = { slope, offset };
+
+        // Monte-Carlo (generates many sets of observations).
+        final int mcRepeat = MONTE_CARLO_RUNS;
+        int mcCount = 0;
+        while (mcCount < mcRepeat) {
+            // Observations.
+            final Point2D.Double[] obs = lineGenerator.generate(numObs);
+
+            final StraightLineProblem problem = new StraightLineProblem(yError);
+            for (int i = 0; i < numObs; i++) {
+                final Point2D.Double p = obs[i];
+                problem.addPoint(p.x, p.y);
+            }
+
+            // Direct solution (using simple regression).
+            final double[] regress = problem.solve();
+
+            // Estimation of the standard deviation (diagonal elements of the
+            // covariance matrix).
+            final PointVectorValuePair optimum
+                = optim.optimize(new MaxEval(Integer.MAX_VALUE),
+                                 problem.getModelFunction(),
+                                 problem.getModelFunctionJacobian(),
+                                 new Target(problem.target()),
+                                 new Weight(problem.weight()),
+                                 new InitialGuess(init));
+            final double[] sigma = optim.computeSigma(optimum.getPoint(), 1e-14);
+
+            // Accumulate statistics.
+            for (int i = 0; i < numParams; i++) {
+                paramsFoundByDirectSolution[i].addValue(regress[i]);
+                sigmaEstimate[i].addValue(sigma[i]);
+            }
+
+            // Next Monte-Carlo.
+            ++mcCount;
+        }
+
+        // Print statistics.
+        final String line = "--------------------------------------------------------------";
+        System.out.println("                 True value       Mean        Std deviation");
+        for (int i = 0; i < numParams; i++) {
+            System.out.println(line);
+            System.out.println("Parameter #" + i);
+
+            StatisticalSummary s = paramsFoundByDirectSolution[i].getSummary();
+            System.out.printf("              %+.6e   %+.6e   %+.6e\n",
+                              init[i],
+                              s.getMean(),
+                              s.getStandardDeviation());
+
+            s = sigmaEstimate[i].getSummary();
+            System.out.printf("sigma: %+.6e (%+.6e)\n",
+                              s.getMean(),
+                              s.getStandardDeviation());
+        }
+        System.out.println(line);
+
+        // Check the error estimation.
+        for (int i = 0; i < numParams; i++) {
+            Assert.assertEquals(paramsFoundByDirectSolution[i].getSummary().getStandardDeviation(),
+                                sigmaEstimate[i].getSummary().getMean(),
+                                8e-2);
+        }
+    }
+
+    /**
+     * In this test, the set of observations is fixed.
+     * Using a Monte-Carlo procedure, it generates sets of parameters,
+     * and determine the parameter change that will result in the
+     * normalized chi-square becoming larger by one than the value from
+     * the best fit solution.
+     * <br/>
+     * The optimization problem solved is defined in class
+     * {@link StraightLineProblem}.
+     * <br/>
+     * The output (on stdout) will be a list of lines containing:
+     * <ul>
+     *  <li>slope of the straight line,</li>
+     *  <li>intercept of the straight line,</li>
+     *  <li>chi-square of the solution defined by the above two values.</li>
+     * </ul>
+     * The output is separated into two blocks (with a blank line between
+     * them); the first block will contain all parameter sets for which
+     * {@code chi2 < chi2_b + 1}
+     * and the second block, all sets for which
+     * {@code chi2 >= chi2_b + 1}
+     * where {@code chi2_b} is the lowest chi-square (corresponding to the
+     * best solution).
+     */
+    @Test
+    public void testParametersErrorMonteCarloParameters() {
+        // Error on the observations.
+        final double yError = 15;
+
+        // True values of the parameters.
+        final double slope = 123.456;
+        final double offset = -98.765;
+
+        // Samples generator.
+        final RandomStraightLinePointGenerator lineGenerator
+            = new RandomStraightLinePointGenerator(slope, offset,
+                                                   yError,
+                                                   -1e3, 1e4,
+                                                   13839013L);
+
+        // Number of observations.
+        final int numObs = 10;
+        // number of parameters.
+        final int numParams = 2;
+
+        // Create a single set of observations.
+        final Point2D.Double[] obs = lineGenerator.generate(numObs);
+
+        final StraightLineProblem problem = new StraightLineProblem(yError);
+        for (int i = 0; i < numObs; i++) {
+            final Point2D.Double p = obs[i];
+            problem.addPoint(p.x, p.y);
+        }
+
+        // Direct solution (using simple regression).
+        final double[] regress = problem.solve();
+
+        // Dummy optimizer (to compute the chi-square).
+        final AbstractLeastSquaresOptimizer optim = new DummyOptimizer();
+        final double[] init = { slope, offset };
+        // Get chi-square of the best parameters set for the given set of
+        // observations.
+        final double bestChi2N = getChi2N(optim, problem, regress);
+        final double[] sigma = optim.computeSigma(regress, 1e-14);
+
+        // Monte-Carlo (generates a grid of parameters).
+        final int mcRepeat = MONTE_CARLO_RUNS;
+        final int gridSize = (int) FastMath.sqrt(mcRepeat);
+
+        // Parameters found for each of Monte-Carlo run.
+        // Index 0 = slope
+        // Index 1 = offset
+        // Index 2 = normalized chi2
+        final List<double[]> paramsAndChi2 = new ArrayList<double[]>(gridSize * gridSize);
+
+        final double slopeRange = 10 * sigma[0];
+        final double offsetRange = 10 * sigma[1];
+        final double minSlope = slope - 0.5 * slopeRange;
+        final double minOffset = offset - 0.5 * offsetRange;
+        final double deltaSlope =  slopeRange/ gridSize;
+        final double deltaOffset = offsetRange / gridSize;
+        for (int i = 0; i < gridSize; i++) {
+            final double s = minSlope + i * deltaSlope;
+            for (int j = 0; j < gridSize; j++) {
+                final double o = minOffset + j * deltaOffset;
+                final double chi2N = getChi2N(optim, problem, new double[] {s, o});
+
+                paramsAndChi2.add(new double[] {s, o, chi2N});
+            }
+        }
+
+        // Output (for use with "gnuplot").
+
+        // Some info.
+
+        // For plotting separately sets of parameters that have a large chi2.
+        final double chi2NPlusOne = bestChi2N + 1;
+        int numLarger = 0;
+
+        final String lineFmt = "%+.10e %+.10e   %.8e\n";
+
+        // Point with smallest chi-square.
+        System.out.printf(lineFmt, regress[0], regress[1], bestChi2N);
+        System.out.println(); // Empty line.
+
+        // Points within the confidence interval.
+        for (double[] d : paramsAndChi2) {
+            if (d[2] <= chi2NPlusOne) {
+                System.out.printf(lineFmt, d[0], d[1], d[2]);
+            }
+        }
+        System.out.println(); // Empty line.
+
+        // Points outside the confidence interval.
+        for (double[] d : paramsAndChi2) {
+            if (d[2] > chi2NPlusOne) {
+                ++numLarger;
+                System.out.printf(lineFmt, d[0], d[1], d[2]);
+            }
+        }
+        System.out.println(); // Empty line.
+
+        System.out.println("# sigma=" + Arrays.toString(sigma));
+        System.out.println("# " + numLarger + " sets filtered out");
+    }
+
+    /**
+     * @return the normalized chi-square.
+     */
+    private double getChi2N(AbstractLeastSquaresOptimizer optim,
+                            StraightLineProblem problem,
+                            double[] params) {
+        final double[] t = problem.target();
+        final double[] w = problem.weight();
+
+        optim.optimize(new MaxEval(Integer.MAX_VALUE),
+                       problem.getModelFunction(),
+                       problem.getModelFunctionJacobian(),
+                       new Target(t),
+                       new Weight(w),
+                       new InitialGuess(params));
+
+        return optim.getChiSquare() / (t.length - params.length);
+    }
+}
+
+/**
+ * A dummy optimizer.
+ * Used for computing the covariance matrix.
+ */
+class DummyOptimizer extends AbstractLeastSquaresOptimizer {
+    public DummyOptimizer() {
+        super(null);
+    }
+
+    /**
+     * This method does nothing and returns a dummy value.
+     */
+    @Override
+    public PointVectorValuePair doOptimize() {
+        final double[] params = getStartPoint();
+        final double[] res = computeResiduals(computeObjectiveValue(params));
+        setCost(computeCost(res));
+        return new PointVectorValuePair(params, null);
+    }
+}

Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optim/nonlinear/vector/jacobian/AbstractLeastSquaresOptimizerTestValidation.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optim/nonlinear/vector/jacobian/CircleProblem.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optim/nonlinear/vector/jacobian/CircleProblem.java?rev=1420684&view=auto
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optim/nonlinear/vector/jacobian/CircleProblem.java (added)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optim/nonlinear/vector/jacobian/CircleProblem.java Wed Dec 12 14:10:38 2012
@@ -0,0 +1,177 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math3.optim.nonlinear.vector.jacobian;
+
+import java.util.ArrayList;
+import org.apache.commons.math3.analysis.MultivariateVectorFunction;
+import org.apache.commons.math3.analysis.MultivariateMatrixFunction;
+import org.apache.commons.math3.util.MathUtils;
+import org.apache.commons.math3.util.FastMath;
+import org.apache.commons.math3.optim.nonlinear.vector.ModelFunction;
+import org.apache.commons.math3.optim.nonlinear.vector.ModelFunctionJacobian;
+
+/**
+ * Class that models a circle.
+ * The parameters of problem are:
+ * <ul>
+ *  <li>the x-coordinate of the circle center,</li>
+ *  <li>the y-coordinate of the circle center,</li>
+ *  <li>the radius of the circle.</li>
+ * </ul>
+ * The model functions are:
+ * <ul>
+ *  <li>for each triplet (cx, cy, r), the (x, y) coordinates of a point on the
+ *   corresponding circle.</li>
+ * </ul>
+ */
+class CircleProblem {
+    /** Cloud of points assumed to be fitted by a circle. */
+    private final ArrayList<double[]> points;
+    /** Error on the x-coordinate of the points. */
+    private final double xSigma;
+    /** Error on the y-coordinate of the points. */
+    private final double ySigma;
+    /** Number of points on the circumference (when searching which
+        model point is closest to a given "observation". */
+    private final int resolution;
+
+    /**
+     * @param xError Assumed error for the x-coordinate of the circle points.
+     * @param yError Assumed error for the y-coordinate of the circle points.
+     * @param searchResolution Number of points to try when searching the one
+     * that is closest to a given "observed" point.
+     */
+    public CircleProblem(double xError,
+                         double yError,
+                         int searchResolution) {
+        points = new ArrayList<double[]>();
+        xSigma = xError;
+        ySigma = yError;
+        resolution = searchResolution;
+    }
+
+    /**
+     * @param xError Assumed error for the x-coordinate of the circle points.
+     * @param yError Assumed error for the y-coordinate of the circle points.
+     */
+    public CircleProblem(double xError,
+                         double yError) {
+        this(xError, yError, 500);
+    }
+
+    public void addPoint(double px, double py) {
+        points.add(new double[] { px, py });
+    }
+
+    public double[] target() {
+        final double[] t = new double[points.size() * 2];
+        for (int i = 0; i < points.size(); i++) {
+            final double[] p = points.get(i);
+            final int index = i * 2;
+            t[index] = p[0];
+            t[index + 1] = p[1];
+        }
+
+        return t;
+    }
+
+    public double[] weight() {
+        final double wX = 1 / (xSigma * xSigma);
+        final double wY = 1 / (ySigma * ySigma);
+        final double[] w = new double[points.size() * 2];
+        for (int i = 0; i < points.size(); i++) {
+            final int index = i * 2;
+            w[index] = wX;
+            w[index + 1] = wY;
+        }
+
+        return w;
+    }
+
+    public ModelFunction getModelFunction() {
+        return new ModelFunction(new MultivariateVectorFunction() {
+                public double[] value(double[] params) {
+                    final double cx = params[0];
+                    final double cy = params[1];
+                    final double r = params[2];
+
+                    final double[] model = new double[points.size() * 2];
+
+                    final double deltaTheta = MathUtils.TWO_PI / resolution;
+                    for (int i = 0; i < points.size(); i++) {
+                        final double[] p = points.get(i);
+                        final double px = p[0];
+                        final double py = p[1];
+
+                        double bestX = 0;
+                        double bestY = 0;
+                        double dMin = Double.POSITIVE_INFINITY;
+
+                        // Find the angle for which the circle passes closest to the
+                        // current point (using a resolution of 100 points along the
+                        // circumference).
+                        for (double theta = 0; theta <= MathUtils.TWO_PI; theta += deltaTheta) {
+                            final double currentX = cx + r * FastMath.cos(theta);
+                            final double currentY = cy + r * FastMath.sin(theta);
+                            final double dX = currentX - px;
+                            final double dY = currentY - py;
+                            final double d = dX * dX + dY * dY;
+                            if (d < dMin) {
+                                dMin = d;
+                                bestX = currentX;
+                                bestY = currentY;
+                            }
+                        }
+
+                        final int index = i * 2;
+                        model[index] = bestX;
+                        model[index + 1] = bestY;
+                    }
+
+                    return model;
+                }
+            });
+    }
+
+    public ModelFunctionJacobian getModelFunctionJacobian() {
+        return new ModelFunctionJacobian(new MultivariateMatrixFunction() {
+                public double[][] value(double[] point) {
+                    return jacobian(point);
+                }
+        });
+    }
+
+    private double[][] jacobian(double[] params) {
+        final double[][] jacobian = new double[points.size() * 2][3];
+
+        for (int i = 0; i < points.size(); i++) {
+            final int index = i * 2;
+            // Partial derivative wrt x-coordinate of center. 
+            jacobian[index][0] = 1;
+            jacobian[index + 1][0] = 0;
+            // Partial derivative wrt y-coordinate of center.
+            jacobian[index][1] = 0;
+            jacobian[index + 1][1] = 1;
+            // Partial derivative wrt radius.
+            final double[] p = points.get(i);
+            jacobian[index][2] = (p[0] - params[0]) / params[2];
+            jacobian[index + 1][2] = (p[1] - params[1]) / params[2];
+        }
+
+        return jacobian;
+    }
+}

Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optim/nonlinear/vector/jacobian/CircleProblem.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optim/nonlinear/vector/jacobian/CircleVectorial.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optim/nonlinear/vector/jacobian/CircleVectorial.java?rev=1420684&view=auto
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optim/nonlinear/vector/jacobian/CircleVectorial.java (added)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optim/nonlinear/vector/jacobian/CircleVectorial.java Wed Dec 12 14:10:38 2012
@@ -0,0 +1,97 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math3.optim.nonlinear.vector.jacobian;
+
+import java.util.ArrayList;
+import org.apache.commons.math3.analysis.MultivariateVectorFunction;
+import org.apache.commons.math3.analysis.MultivariateMatrixFunction;
+import org.apache.commons.math3.geometry.euclidean.twod.Vector2D;
+import org.apache.commons.math3.optim.nonlinear.vector.ModelFunction;
+import org.apache.commons.math3.optim.nonlinear.vector.ModelFunctionJacobian;
+
+/**
+ * Class used in the tests.
+ */
+class CircleVectorial {
+    private ArrayList<Vector2D> points;
+
+    public CircleVectorial() {
+        points  = new ArrayList<Vector2D>();
+    }
+
+    public void addPoint(double px, double py) {
+        points.add(new Vector2D(px, py));
+    }
+
+    public int getN() {
+        return points.size();
+    }
+
+    public double getRadius(Vector2D center) {
+        double r = 0;
+        for (Vector2D point : points) {
+            r += point.distance(center);
+        }
+        return r / points.size();
+    }
+
+    public ModelFunction getModelFunction() {
+        return new ModelFunction(new MultivariateVectorFunction() {
+                public double[] value(double[] params) {
+                    Vector2D center = new Vector2D(params[0], params[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;
+                }
+        });
+    }
+
+    public ModelFunctionJacobian getModelFunctionJacobian() {
+        return new ModelFunctionJacobian(new MultivariateMatrixFunction() {
+                public double[][] value(double[] params) {
+                    final int n = points.size();
+                    final Vector2D center = new Vector2D(params[0], params[1]);
+
+                    double dRdX = 0;
+                    double dRdY = 0;
+                    for (Vector2D pk : points) {
+                        double dk = pk.distance(center);
+                        dRdX += (center.getX() - pk.getX()) / dk;
+                        dRdY += (center.getY() - pk.getY()) / dk;
+                    }
+                    dRdX /= n;
+                    dRdY /= n;
+
+                    // Jacobian of the radius residuals.
+                    double[][] jacobian = new double[n][2];
+                    for (int i = 0; i < n; i++) {
+                        final Vector2D pi = points.get(i);
+                        final double di = pi.distance(center);
+                        jacobian[i][0] = (center.getX() - pi.getX()) / di - dRdX;
+                        jacobian[i][1] = (center.getY() - pi.getY()) / di - dRdY;
+                    }
+
+                    return jacobian;
+                }
+        });
+    }
+}

Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optim/nonlinear/vector/jacobian/CircleVectorial.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optim/nonlinear/vector/jacobian/GaussNewtonOptimizerTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optim/nonlinear/vector/jacobian/GaussNewtonOptimizerTest.java?rev=1420684&view=auto
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optim/nonlinear/vector/jacobian/GaussNewtonOptimizerTest.java (added)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optim/nonlinear/vector/jacobian/GaussNewtonOptimizerTest.java Wed Dec 12 14:10:38 2012
@@ -0,0 +1,159 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math3.optim.nonlinear.vector.jacobian;
+
+import java.io.IOException;
+import org.apache.commons.math3.exception.ConvergenceException;
+import org.apache.commons.math3.exception.TooManyEvaluationsException;
+import org.apache.commons.math3.optim.SimpleVectorValueChecker;
+import org.apache.commons.math3.optim.InitialGuess;
+import org.apache.commons.math3.optim.MaxEval;
+import org.apache.commons.math3.optim.nonlinear.vector.Target;
+import org.apache.commons.math3.optim.nonlinear.vector.Weight;
+import org.apache.commons.math3.optim.nonlinear.vector.ModelFunction;
+import org.apache.commons.math3.optim.nonlinear.vector.ModelFunctionJacobian;
+import org.junit.Test;
+
+/**
+ * <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
+ * href="http://www.netlib.org/minpack/ex/file22">file22</a> test files.
+ * The redistribution policy for MINPACK is available <a
+ * href="http://www.netlib.org/minpack/disclaimer">here</a>, for
+ * convenience, it is reproduced below.</p>
+
+ * <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0">
+ * <tr><td>
+ *    Minpack Copyright Notice (1999) University of Chicago.
+ *    All rights reserved
+ * </td></tr>
+ * <tr><td>
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * <ol>
+ *  <li>Redistributions of source code must retain the above copyright
+ *      notice, this list of conditions and the following disclaimer.</li>
+ * <li>Redistributions in binary form must reproduce the above
+ *     copyright notice, this list of conditions and the following
+ *     disclaimer in the documentation and/or other materials provided
+ *     with the distribution.</li>
+ * <li>The end-user documentation included with the redistribution, if any,
+ *     must include the following acknowledgment:
+ *     <code>This product includes software developed by the University of
+ *           Chicago, as Operator of Argonne National Laboratory.</code>
+ *     Alternately, this acknowledgment may appear in the software itself,
+ *     if and wherever such third-party acknowledgments normally appear.</li>
+ * <li><strong>WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS"
+ *     WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE
+ *     UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND
+ *     THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR
+ *     IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES
+ *     OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE
+ *     OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY
+ *     OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR
+ *     USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF
+ *     THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4)
+ *     DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION
+ *     UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL
+ *     BE CORRECTED.</strong></li>
+ * <li><strong>LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT
+ *     HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF
+ *     ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT,
+ *     INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF
+ *     ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF
+ *     PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER
+ *     SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT
+ *     (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE,
+ *     EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE
+ *     POSSIBILITY OF SUCH LOSS OR DAMAGES.</strong></li>
+ * <ol></td></tr>
+ * </table>
+
+ * @author Argonne National Laboratory. MINPACK project. March 1980 (original fortran minpack tests)
+ * @author Burton S. Garbow (original fortran minpack tests)
+ * @author Kenneth E. Hillstrom (original fortran minpack tests)
+ * @author Jorge J. More (original fortran minpack tests)
+ * @author Luc Maisonobe (non-minpack tests and minpack tests Java translation)
+ */
+public class GaussNewtonOptimizerTest
+    extends AbstractLeastSquaresOptimizerAbstractTest {
+
+    @Override
+    public AbstractLeastSquaresOptimizer createOptimizer() {
+        return new GaussNewtonOptimizer(new SimpleVectorValueChecker(1.0e-6, 1.0e-6));
+    }
+
+    @Override
+    @Test(expected = ConvergenceException.class)
+    public void testMoreEstimatedParametersSimple() {
+        /*
+         * Exception is expected with this optimizer
+         */
+        super.testMoreEstimatedParametersSimple();
+    }
+
+    @Override
+    @Test(expected=ConvergenceException.class)
+    public void testMoreEstimatedParametersUnsorted() {
+        /*
+         * Exception is expected with this optimizer
+         */
+        super.testMoreEstimatedParametersUnsorted();
+    }
+
+    @Test(expected=TooManyEvaluationsException.class)
+    public void testMaxEvaluations() throws Exception {
+        CircleVectorial circle = new CircleVectorial();
+        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);
+
+        GaussNewtonOptimizer optimizer
+            = new GaussNewtonOptimizer(new SimpleVectorValueChecker(1e-30, 1e-30));
+
+        optimizer.optimize(new MaxEval(100),
+                           circle.getModelFunction(),
+                           circle.getModelFunctionJacobian(),
+                           new Target(new double[] { 0, 0, 0, 0, 0 }),
+                           new Weight(new double[] { 1, 1, 1, 1, 1 }),
+                           new InitialGuess(new double[] { 98.680, 47.345 }));
+    }
+
+    @Override
+    @Test(expected=ConvergenceException.class)
+    public void testCircleFittingBadInit() {
+        /*
+         * This test does not converge with this optimizer.
+         */
+        super.testCircleFittingBadInit();
+    }
+
+    @Override
+    @Test(expected = ConvergenceException.class)
+    public void testHahn1()
+        throws IOException {
+        /*
+         * TODO This test leads to a singular problem with the Gauss-Newton
+         * optimizer. This should be inquired.
+         */
+        super.testHahn1();
+    }
+}

Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optim/nonlinear/vector/jacobian/GaussNewtonOptimizerTest.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optim/nonlinear/vector/jacobian/LevenbergMarquardtOptimizerTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optim/nonlinear/vector/jacobian/LevenbergMarquardtOptimizerTest.java?rev=1420684&view=auto
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optim/nonlinear/vector/jacobian/LevenbergMarquardtOptimizerTest.java (added)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optim/nonlinear/vector/jacobian/LevenbergMarquardtOptimizerTest.java Wed Dec 12 14:10:38 2012
@@ -0,0 +1,403 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math3.optim.nonlinear.vector.jacobian;
+
+import java.io.Serializable;
+import java.util.ArrayList;
+import java.util.List;
+import org.apache.commons.math3.optim.PointVectorValuePair;
+import org.apache.commons.math3.optim.InitialGuess;
+import org.apache.commons.math3.optim.MaxEval;
+import org.apache.commons.math3.optim.nonlinear.vector.Target;
+import org.apache.commons.math3.optim.nonlinear.vector.Weight;
+import org.apache.commons.math3.optim.nonlinear.vector.ModelFunction;
+import org.apache.commons.math3.optim.nonlinear.vector.ModelFunctionJacobian;
+import org.apache.commons.math3.analysis.MultivariateVectorFunction;
+import org.apache.commons.math3.analysis.MultivariateMatrixFunction;
+import org.apache.commons.math3.exception.ConvergenceException;
+import org.apache.commons.math3.exception.DimensionMismatchException;
+import org.apache.commons.math3.exception.TooManyEvaluationsException;
+import org.apache.commons.math3.geometry.euclidean.twod.Vector2D;
+import org.apache.commons.math3.linear.SingularMatrixException;
+import org.apache.commons.math3.util.FastMath;
+import org.apache.commons.math3.util.Precision;
+import org.junit.Assert;
+import org.junit.Test;
+import org.junit.Ignore;
+
+/**
+ * <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
+ * href="http://www.netlib.org/minpack/ex/file22">file22</a> test files.
+ * The redistribution policy for MINPACK is available <a
+ * href="http://www.netlib.org/minpack/disclaimer">here</a>, for
+ * convenience, it is reproduced below.</p>
+
+ * <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0">
+ * <tr><td>
+ *    Minpack Copyright Notice (1999) University of Chicago.
+ *    All rights reserved
+ * </td></tr>
+ * <tr><td>
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * <ol>
+ *  <li>Redistributions of source code must retain the above copyright
+ *      notice, this list of conditions and the following disclaimer.</li>
+ * <li>Redistributions in binary form must reproduce the above
+ *     copyright notice, this list of conditions and the following
+ *     disclaimer in the documentation and/or other materials provided
+ *     with the distribution.</li>
+ * <li>The end-user documentation included with the redistribution, if any,
+ *     must include the following acknowledgment:
+ *     <code>This product includes software developed by the University of
+ *           Chicago, as Operator of Argonne National Laboratory.</code>
+ *     Alternately, this acknowledgment may appear in the software itself,
+ *     if and wherever such third-party acknowledgments normally appear.</li>
+ * <li><strong>WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS"
+ *     WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE
+ *     UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND
+ *     THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR
+ *     IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES
+ *     OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE
+ *     OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY
+ *     OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR
+ *     USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF
+ *     THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4)
+ *     DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION
+ *     UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL
+ *     BE CORRECTED.</strong></li>
+ * <li><strong>LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT
+ *     HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF
+ *     ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT,
+ *     INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF
+ *     ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF
+ *     PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER
+ *     SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT
+ *     (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE,
+ *     EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE
+ *     POSSIBILITY OF SUCH LOSS OR DAMAGES.</strong></li>
+ * <ol></td></tr>
+ * </table>
+
+ * @author Argonne National Laboratory. MINPACK project. March 1980 (original fortran minpack tests)
+ * @author Burton S. Garbow (original fortran minpack tests)
+ * @author Kenneth E. Hillstrom (original fortran minpack tests)
+ * @author Jorge J. More (original fortran minpack tests)
+ * @author Luc Maisonobe (non-minpack tests and minpack tests Java translation)
+ */
+public class LevenbergMarquardtOptimizerTest
+    extends AbstractLeastSquaresOptimizerAbstractTest {
+    @Override
+    public AbstractLeastSquaresOptimizer createOptimizer() {
+        return new LevenbergMarquardtOptimizer();
+    }
+
+    @Override
+    @Test(expected=SingularMatrixException.class)
+    public void testNonInvertible() {
+        /*
+         * Overrides the method from parent class, since the default singularity
+         * threshold (1e-14) does not trigger the expected exception.
+         */
+        LinearProblem problem = new LinearProblem(new double[][] {
+                {  1, 2, -3 },
+                {  2, 1,  3 },
+                { -3, 0, -9 }
+        }, new double[] { 1, 1, 1 });
+
+        AbstractLeastSquaresOptimizer optimizer = createOptimizer();
+        PointVectorValuePair optimum
+            = optimizer.optimize(new MaxEval(100),
+                                 problem.getModelFunction(),
+                                 problem.getModelFunctionJacobian(),
+                                 problem.getTarget(),
+                                 new Weight(new double[] { 1, 1, 1 }),
+                                 new InitialGuess(new double[] { 0, 0, 0 }));
+        Assert.assertTrue(FastMath.sqrt(optimizer.getTargetSize()) * optimizer.getRMS() > 0.6);
+
+        optimizer.computeCovariances(optimum.getPoint(), 1.5e-14);
+    }
+
+    @Test
+    public void testControlParameters() {
+        CircleVectorial circle = new CircleVectorial();
+        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.getModelFunction(),
+                      circle.getModelFunctionJacobian(),
+                      0.1, 10, 1.0e-14, 1.0e-16, 1.0e-10, false);
+        checkEstimate(circle.getModelFunction(),
+                      circle.getModelFunctionJacobian(),
+                      0.1, 10, 1.0e-15, 1.0e-17, 1.0e-10, true);
+        checkEstimate(circle.getModelFunction(),
+                      circle.getModelFunctionJacobian(),
+                      0.1,  5, 1.0e-15, 1.0e-16, 1.0e-10, true);
+        circle.addPoint(300, -300);
+        checkEstimate(circle.getModelFunction(),
+                      circle.getModelFunctionJacobian(),
+                      0.1, 20, 1.0e-18, 1.0e-16, 1.0e-10, true);
+    }
+
+    private void checkEstimate(ModelFunction problem,
+                               ModelFunctionJacobian problemJacobian,
+                               double initialStepBoundFactor, int maxCostEval,
+                               double costRelativeTolerance, double parRelativeTolerance,
+                               double orthoTolerance, boolean shouldFail) {
+        try {
+            LevenbergMarquardtOptimizer optimizer
+                = new LevenbergMarquardtOptimizer(initialStepBoundFactor,
+                                                  costRelativeTolerance,
+                                                  parRelativeTolerance,
+                                                  orthoTolerance,
+                                                  Precision.SAFE_MIN);
+            optimizer.optimize(new MaxEval(maxCostEval),
+                               problem,
+                               problemJacobian,
+                               new Target(new double[] { 0, 0, 0, 0, 0 }),
+                               new Weight(new double[] { 1, 1, 1, 1, 1 }),
+                               new InitialGuess(new double[] { 98.680, 47.345 }));
+            Assert.assertTrue(!shouldFail);
+        } catch (DimensionMismatchException ee) {
+            Assert.assertTrue(shouldFail);
+        } catch (TooManyEvaluationsException ee) {
+            Assert.assertTrue(shouldFail);
+        }
+    }
+
+    /**
+     * Non-linear test case: fitting of decay curve (from Chapter 8 of
+     * Bevington's textbook, "Data reduction and analysis for the physical sciences").
+     * XXX The expected ("reference") values may not be accurate and the tolerance too
+     * relaxed for this test to be currently really useful (the issue is under
+     * investigation).
+     */
+    @Test
+    public void testBevington() {
+        final double[][] dataPoints = {
+            // column 1 = times
+            { 15, 30, 45, 60, 75, 90, 105, 120, 135, 150,
+              165, 180, 195, 210, 225, 240, 255, 270, 285, 300,
+              315, 330, 345, 360, 375, 390, 405, 420, 435, 450,
+              465, 480, 495, 510, 525, 540, 555, 570, 585, 600,
+              615, 630, 645, 660, 675, 690, 705, 720, 735, 750,
+              765, 780, 795, 810, 825, 840, 855, 870, 885, },
+            // column 2 = measured counts
+            { 775, 479, 380, 302, 185, 157, 137, 119, 110, 89,
+              74, 61, 66, 68, 48, 54, 51, 46, 55, 29,
+              28, 37, 49, 26, 35, 29, 31, 24, 25, 35,
+              24, 30, 26, 28, 21, 18, 20, 27, 17, 17,
+              14, 17, 24, 11, 22, 17, 12, 10, 13, 16,
+              9, 9, 14, 21, 17, 13, 12, 18, 10, },
+        };
+
+        final BevingtonProblem problem = new BevingtonProblem();
+
+        final int len = dataPoints[0].length;
+        final double[] weights = new double[len];
+        for (int i = 0; i < len; i++) {
+            problem.addPoint(dataPoints[0][i],
+                             dataPoints[1][i]);
+
+            weights[i] = 1 / dataPoints[1][i];
+        }
+
+        final LevenbergMarquardtOptimizer optimizer
+            = new LevenbergMarquardtOptimizer();
+
+        final PointVectorValuePair optimum
+            = optimizer.optimize(new MaxEval(100),
+                                 problem.getModelFunction(),
+                                 problem.getModelFunctionJacobian(),
+                                 new Target(dataPoints[1]),
+                                 new Weight(weights),
+                                 new InitialGuess(new double[] { 10, 900, 80, 27, 225 }));
+
+        final double[] solution = optimum.getPoint();
+        final double[] expectedSolution = { 10.4, 958.3, 131.4, 33.9, 205.0 };
+
+        final double[][] covarMatrix = optimizer.computeCovariances(solution, 1e-14);
+        final double[][] expectedCovarMatrix = {
+            { 3.38, -3.69, 27.98, -2.34, -49.24 },
+            { -3.69, 2492.26, 81.89, -69.21, -8.9 },
+            { 27.98, 81.89, 468.99, -44.22, -615.44 },
+            { -2.34, -69.21, -44.22, 6.39, 53.80 },
+            { -49.24, -8.9, -615.44, 53.8, 929.45 }
+        };
+
+        final int numParams = expectedSolution.length;
+
+        // Check that the computed solution is within the reference error range.
+        for (int i = 0; i < numParams; i++) {
+            final double error = FastMath.sqrt(expectedCovarMatrix[i][i]);
+            Assert.assertEquals("Parameter " + i, expectedSolution[i], solution[i], error);
+        }
+
+        // Check that each entry of the computed covariance matrix is within 10%
+        // of the reference matrix entry.
+        for (int i = 0; i < numParams; i++) {
+            for (int j = 0; j < numParams; j++) {
+                Assert.assertEquals("Covariance matrix [" + i + "][" + j + "]",
+                                    expectedCovarMatrix[i][j],
+                                    covarMatrix[i][j],
+                                    FastMath.abs(0.1 * expectedCovarMatrix[i][j]));
+            }
+        }
+    }
+
+    @Test
+    public void testCircleFitting2() {
+        final double xCenter = 123.456;
+        final double yCenter = 654.321;
+        final double xSigma = 10;
+        final double ySigma = 15;
+        final double radius = 111.111;
+        // The test is extremely sensitive to the seed.
+        final long seed = 59421061L;
+        final RandomCirclePointGenerator factory
+            = new RandomCirclePointGenerator(xCenter, yCenter, radius,
+                                             xSigma, ySigma,
+                                             seed);
+        final CircleProblem circle = new CircleProblem(xSigma, ySigma);
+
+        final int numPoints = 10;
+        for (Vector2D p : factory.generate(numPoints)) {
+            circle.addPoint(p.getX(), p.getY());
+        }
+
+        // First guess for the center's coordinates and radius.
+        final double[] init = { 90, 659, 115 };
+
+        final LevenbergMarquardtOptimizer optimizer
+            = new LevenbergMarquardtOptimizer();
+        final PointVectorValuePair optimum = optimizer.optimize(new MaxEval(100),
+                                                                circle.getModelFunction(),
+                                                                circle.getModelFunctionJacobian(),
+                                                                new Target(circle.target()),
+                                                                new Weight(circle.weight()),
+                                                                new InitialGuess(init));
+
+        final double[] paramFound = optimum.getPoint();
+
+        // Retrieve errors estimation.
+        final double[] asymptoticStandardErrorFound = optimizer.computeSigma(paramFound, 1e-14);
+
+        // Check that the parameters are found within the assumed error bars.
+        Assert.assertEquals(xCenter, paramFound[0], asymptoticStandardErrorFound[0]);
+        Assert.assertEquals(yCenter, paramFound[1], asymptoticStandardErrorFound[1]);
+        Assert.assertEquals(radius, paramFound[2], asymptoticStandardErrorFound[2]);
+    }
+
+    private static class QuadraticProblem {
+        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 ModelFunction getModelFunction() {
+            return new ModelFunction(new MultivariateVectorFunction() {
+                    public double[] value(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 ModelFunctionJacobian getModelFunctionJacobian() {
+            return new ModelFunctionJacobian(new MultivariateMatrixFunction() {
+                    public double[][] value(double[] params) {                    
+                        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;
+                    }
+                });
+        }
+    }
+
+    private static class BevingtonProblem {
+        private List<Double> time;
+        private List<Double> count;
+
+        public BevingtonProblem() {
+            time = new ArrayList<Double>();
+            count = new ArrayList<Double>();
+        }
+
+        public void addPoint(double t, double c) {
+            time.add(t);
+            count.add(c);
+        }
+
+        public ModelFunction getModelFunction() {
+            return new ModelFunction(new MultivariateVectorFunction() {
+                    public double[] value(double[] params) {
+                        double[] values = new double[time.size()];
+                        for (int i = 0; i < values.length; ++i) {
+                            final double t = time.get(i);
+                            values[i] = params[0] +
+                                params[1] * Math.exp(-t / params[3]) +
+                                params[2] * Math.exp(-t / params[4]);
+                        }
+                        return values;
+                    }
+                });
+        }
+
+        public ModelFunctionJacobian getModelFunctionJacobian() {
+            return new ModelFunctionJacobian(new MultivariateMatrixFunction() {
+                    public double[][] value(double[] params) {
+                        double[][] jacobian = new double[time.size()][5];
+
+                        for (int i = 0; i < jacobian.length; ++i) {
+                            final double t = time.get(i);
+                            jacobian[i][0] = 1;
+
+                            final double p3 =  params[3];
+                            final double p4 =  params[4];
+                            final double tOp3 = t / p3;
+                            final double tOp4 = t / p4;
+                            jacobian[i][1] = Math.exp(-tOp3);
+                            jacobian[i][2] = Math.exp(-tOp4);
+                            jacobian[i][3] = params[1] * Math.exp(-tOp3) * tOp3 / p3;
+                            jacobian[i][4] = params[2] * Math.exp(-tOp4) * tOp4 / p4;
+                        }
+                        return jacobian;
+                    }
+                });
+        }
+    }
+}

Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optim/nonlinear/vector/jacobian/LevenbergMarquardtOptimizerTest.java
------------------------------------------------------------------------------
    svn:eol-style = native



Mime
View raw message