aboutsummaryrefslogtreecommitdiff
path: root/contrib/cone-vector.py
diff options
context:
space:
mode:
Diffstat (limited to 'contrib/cone-vector.py')
-rw-r--r--contrib/cone-vector.py24
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