You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by er...@apache.org on 2012/12/12 15:11:04 UTC
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...
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