Implement rejection sampling based cone sampling.

* contrib/cone-vector.py: Import log.
(random_planar_angle_pdf, random_vector_on_spherical_cap_pdf): New
functions.

1 files changed, 23 insertions, 1 deletions

diff --git a/contrib/cone-vector.py b/contrib/cone-vector.py
index 65e08ee..1f1fb1d 100644
--- a/contrib/cone-vector.py
+++ b/contrib/cone-vector.py
@@ -16,7 +16,7 @@
 # along with this program.  If not, see
 # <https://www.gnu.org/licenses/>.
 
-from numpy import arcsin, cos, dot, empty, ones, sin, sqrt, pi, where
+from numpy import arcsin, cos, dot, empty, log, ones, sin, sqrt, pi, where
 from numpy.random import randn, random
 from numpy.linalg import norm
 from scipy.special import betainc, betaincinv
@@ -62,6 +62,17 @@ def random_planar_angle_cdf(maximum_planar_angle, dim):
         planar_angle2solid_angle_fraction(maximum_planar_angle, dim)*random(),
         dim)
 
+def random_planar_angle_pdf(maximum_planar_angle, dim):
+    """Return a random planar angle using rejection sampling."""
+    # We apply the log function just to prevent the floats from
+    # underflowing.
+    box_height = (dim-2)*log(sin(min(maximum_planar_angle, pi/2)))
+    while True:
+        theta = maximum_planar_angle*random()
+        f = box_height + log(random())
+        if f < (dim-2)*log(sin(theta)):
+            return theta
+
 def random_vector_on_disk(axis, planar_angle):
     """Return a random vector uniformly distributed on the periphery of a
 disk."""
@@ -77,6 +88,17 @@ cap. The random planar angle is generated using inverse transform
 sampling."""
     return random_vector_on_disk(axis, random_planar_angle_cdf(maximum_planar_angle, axis.size))
 
+def random_vector_on_spherical_cap_pdf(axis, maximum_planar_angle):
+    """Return a random vector uniformly distributed on a spherical
+cap. The random planar angle is generated using rejection sampling.
+
+This function is more numerically stable than
+random_vector_on_spherical_cap_cdf for large dimensions and small
+angles.
+
+    """
+    return random_vector_on_disk(axis, random_planar_angle_pdf(maximum_planar_angle, axis.size))
+
 def sample_code():
     """Run some sample code testing the defined functions."""
     dim = 100