MATH-931

Efficiency improvement and unit test (thanks to Sean Owen).


git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1440734 13f79535-47bb-0310-9956-ffa450edef68

src/changes/changes.xml

@@ -55,6 +55,9 @@ This is a minor release: It combines bug fixes and new features.
   Changes to existing features were made in a backwards-compatible
   way such as to allow drop-in replacement of the v3.1[.1] JAR file.
 ">
+      <action dev="erans" type="update" issue="MATH-931" due-to="Sean Owen">
+        Greater efficiency in "UnitSphereRandomVectorGenerator".
+      </action>
       <action dev="tn" type="fix" issue="MATH-930">
         Improved class javadoc wrt convergence criteria and added 
         additional constructors to override the default epsilon and cut-off

src/main/java/org/apache/commons/math3/random/UnitSphereRandomVectorGenerator.java

@@ -59,18 +59,17 @@ public class UnitSphereRandomVectorGenerator
 
     /** {@inheritDoc} */
     public double[] nextVector() {
-
         final double[] v = new double[dimension];
 
-        double normSq;
-        do {
-            normSq = 0;
-            for (int i = 0; i < dimension; i++) {
-                final double comp = 2 * rand.nextDouble() - 1;
-                v[i] = comp;
-                normSq += comp * comp;
-            }
-        } while (normSq > 1);
+        // See http://mathworld.wolfram.com/SpherePointPicking.html for example.
+        // Pick a point by choosing a standard Gaussian for each element, and then 
+        // normalizing to unit length.
+        double normSq = 0;
+        for (int i = 0; i < dimension; i++) {
+            final double comp = rand.nextGaussian();
+            v[i] = comp;
+            normSq += comp * comp;
+        }
 
         final double f = 1 / FastMath.sqrt(normSq);
         for (int i = 0; i < dimension; i++) {
@@ -78,7 +77,5 @@ public class UnitSphereRandomVectorGenerator
         }
 
         return v;
-
     }
-
 }

src/test/java/org/apache/commons/math3/random/UnitSphereRandomVectorGeneratorTest.java

@@ -0,0 +1,68 @@
+/*
+ * 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.random;
+
+import org.apache.commons.math3.util.FastMath;
+import org.junit.Assert;
+import org.junit.Test;
+
+public class UnitSphereRandomVectorGeneratorTest {
+    /**
+     * Test the distribution of points from {@link UnitSphereRandomVectorGenerator#nextVector()}
+     * in two dimensions.
+     */
+    @Test
+    public void test2DDistribution() {
+      
+        RandomGenerator rg = new JDKRandomGenerator();
+        rg.setSeed(17399225432l);
+        UnitSphereRandomVectorGenerator generator = new UnitSphereRandomVectorGenerator(2, rg);
+
+        // In 2D, angles with a given vector should be uniformly distributed
+        int[] angleBuckets = new int[100];
+        int steps = 1000000;
+        for (int i = 0; i < steps; ++i) {
+            final double[] v = generator.nextVector();
+            Assert.assertEquals(2, v.length);
+            Assert.assertEquals(1, length(v), 1e-10);
+            // Compute angle formed with vector (1,0)
+            // Cosine of angle is their dot product, because both are unit length
+            // Dot product here is just the first element of the vector by construction
+            final double angle = FastMath.acos(v[0]);
+            final int bucket = (int) (angleBuckets.length * (angle / FastMath.PI));
+            ++angleBuckets[bucket];
+        }
+
+        // Simplistic test for roughly even distribution
+        final int expectedBucketSize = steps / angleBuckets.length;
+        for (int bucket : angleBuckets) {
+            Assert.assertTrue("Bucket count " + bucket + " vs expected " + expectedBucketSize,
+                              FastMath.abs(expectedBucketSize - bucket) < 350);
+        }
+    }
+
+    /**
+     * @return length (L2 norm) of given vector
+     */
+    private static double length(double[] vector) {
+        double total = 0;
+        for (double d : vector) {
+            total += d * d;
+        }
+        return FastMath.sqrt(total);
+    }
+}