diff options
| -rw-r--r-- | contrib/cone-vector.py | 24 |
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 |