commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From er...@apache.org
Subject svn commit: r1330801 - in /commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general: CircleProblem.java LevenbergMarquardtOptimizerTest.java RandomCirclePointGenerator.java
Date Thu, 26 Apr 2012 12:13:49 GMT
Author: erans
Date: Thu Apr 26 12:13:49 2012
New Revision: 1330801

URL: http://svn.apache.org/viewvc?rev=1330801&view=rev
Log:
Added unit test.

Added:
    commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/CircleProblem.java
  (with props)
    commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/RandomCirclePointGenerator.java
  (with props)
Modified:
    commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/LevenbergMarquardtOptimizerTest.java

Added: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/CircleProblem.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/CircleProblem.java?rev=1330801&view=auto
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/CircleProblem.java
(added)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/CircleProblem.java
Thu Apr 26 12:13:49 2012
@@ -0,0 +1,173 @@
+/*
+ * 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.optimization.general;
+
+import java.util.Arrays;
+import java.util.ArrayList;
+import org.apache.commons.math3.analysis.DifferentiableMultivariateVectorFunction;
+import org.apache.commons.math3.analysis.MultivariateMatrixFunction;
+import org.apache.commons.math3.util.MathUtils;
+import org.apache.commons.math3.util.FastMath;
+
+/**
+ * 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 implements DifferentiableMultivariateVectorFunction {
+    /** 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 xError 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 xError 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 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 MultivariateMatrixFunction jacobian() {
+        return 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/optimization/general/CircleProblem.java
------------------------------------------------------------------------------
    svn:eol-style = native

Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/LevenbergMarquardtOptimizerTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/LevenbergMarquardtOptimizerTest.java?rev=1330801&r1=1330800&r2=1330801&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/LevenbergMarquardtOptimizerTest.java
(original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/LevenbergMarquardtOptimizerTest.java
Thu Apr 26 12:13:49 2012
@@ -590,6 +590,55 @@ public class LevenbergMarquardtOptimizer
         }
     }
 
+    @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;
+        final RandomCirclePointGenerator factory
+            = new RandomCirclePointGenerator(xCenter, yCenter, radius,
+                                             xSigma, ySigma,
+                                             59421063L);
+        final CircleProblem circle = new CircleProblem(xSigma, ySigma);
+
+        final int numPoints = 10;
+        for (Point2D.Double p : factory.generate(numPoints)) {
+            circle.addPoint(p.x, p.y);
+            // System.out.println(p.x + " " + p.y);
+        }
+
+        // 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(100, circle,
+                                                                circle.target(), circle.weight(),
+                                                                init);
+
+        final double[] paramFound = optimum.getPoint();
+
+        // Retrieve errors estimation.
+        final double[][] covMatrix = optimizer.getCovariances();
+        final double[] asymptoticStandardErrorFound = optimizer.guessParametersErrors();
+        final double[] sigmaFound = new double[covMatrix.length];
+        for (int i = 0; i < covMatrix.length; i++) {
+            sigmaFound[i] = FastMath.sqrt(covMatrix[i][i]);
+//             System.out.println("i=" + i + " value=" + paramFound[i]
+//                                + " sigma=" + sigmaFound[i]
+//                                + " ase=" + asymptoticStandardErrorFound[i]);
+        }
+
+        // System.out.println("chi2=" + optimizer.getChiSquare());
+
+        // 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 LinearProblem implements DifferentiableMultivariateVectorFunction,
Serializable {
 
         private static final long serialVersionUID = 703247177355019415L;

Added: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/RandomCirclePointGenerator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/RandomCirclePointGenerator.java?rev=1330801&view=auto
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/RandomCirclePointGenerator.java
(added)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/RandomCirclePointGenerator.java
Thu Apr 26 12:13:49 2012
@@ -0,0 +1,95 @@
+/*
+ * 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.optimization.general;
+
+import java.awt.geom.Point2D;
+import org.apache.commons.math3.random.RandomData;
+import org.apache.commons.math3.random.RandomDataImpl;
+import org.apache.commons.math3.random.Well44497b;
+import org.apache.commons.math3.util.MathUtils;
+import org.apache.commons.math3.util.FastMath;
+
+/**
+ * Factory for generating a cloud of points that approximate a circle.
+ */
+public class RandomCirclePointGenerator {
+    /** RNG. */
+    private final RandomData random;
+    /** Radius of the circle. */
+    private final double radius;
+    /** x-coordinate of the circle center. */
+    private final double x;
+    /** y-coordinate of the circle center. */
+    private final double y;
+    /** Error on the x-coordinate of the center. */
+    private final double xSigma;
+    /** Error on the x-coordinate of the center. */
+    private final double ySigma;
+
+    /**
+     * @param x Abscissa of the circle center.
+     * @param y Ordinate of the circle center.
+     * @param radius Radius of the circle.
+     * @param xSigma Error on the x-coordinate of the circumferenc points.
+     * @param ySigma Error on the y-coordinate of the circumferenc points.
+     * @param seed RNG seed.
+     */
+    public RandomCirclePointGenerator(double x,
+                                      double y,
+                                      double radius,
+                                      double xSigma,
+                                      double ySigma,
+                                      long seed) {
+        random = new RandomDataImpl(new Well44497b((seed)));
+        this.radius = radius;
+        this.x = x;
+        this.y = y;
+        this.xSigma = xSigma;
+        this.ySigma = ySigma;
+    }
+
+    /**
+     * Point generator.
+     *
+     * @param n Number of points to create.
+     * @return the cloud of {@code n} points.
+     */
+    public Point2D.Double[] generate(int n) {
+        final Point2D.Double[] cloud = new Point2D.Double[n];
+        for (int i = 0; i < n; i++) {
+            cloud[i] = create();
+        }
+        return cloud;
+    }
+
+    /**
+     * Create one point.
+     *
+     * @return a point.
+     */
+    private Point2D.Double create() {
+        final double cX = random.nextGaussian(x, xSigma);
+        final double cY = random.nextGaussian(y, ySigma);
+        final double t = random.nextUniform(0, MathUtils.TWO_PI);
+
+        final double pX = cX + radius * FastMath.cos(t);
+        final double pY = cY + radius * FastMath.sin(t);
+
+        return new Point2D.Double(pX, pY);
+    }
+}

Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/RandomCirclePointGenerator.java
------------------------------------------------------------------------------
    svn:eol-style = native



Mime
View raw message