From 7fdcd4d355274bdd70fa1de0ff9723790a45634f Mon Sep 17 00:00:00 2001 From: Arun Isaac Date: Fri, 26 Mar 2021 14:13:48 +0530 Subject: Remove inverse transform sampling based functions. * src/sambal/sambal.py: Do not import numpy.arcsin, numpy.ones, numpy.where, scipy.special.betainc, scipy.special.betaincinv. (planar_angle2solid_angle_fraction, solid_angle_fraction2planar_angle, random_planar_angle_cdf, random_vector_on_spherical_cap_cdf): Delete functions. (random_vector_on_spherical_cap_pdf): Remove reference to random_vector_on_spherical_cap_cdf in docstring. Rename to random_on_cap. --- src/sambal/sambal.py | 37 ++----------------------------------- 1 file changed, 2 insertions(+), 35 deletions(-) diff --git a/src/sambal/sambal.py b/src/sambal/sambal.py index de68d1c..3a35438 100644 --- a/src/sambal/sambal.py +++ b/src/sambal/sambal.py @@ -16,32 +16,15 @@ # along with this program. If not, see # . -from numpy import arcsin, cos, dot, empty, log, ones, sin, sqrt, pi, where +from numpy import cos, dot, empty, log, sin, sqrt, pi from numpy.random import randn, random from numpy.linalg import norm -from scipy.special import betainc, betaincinv def random_on_sphere(dim): """Return a random vector uniformly distributed on the unit sphere.""" x = randn(dim) return x / norm(x) -def planar_angle2solid_angle_fraction(planar_angle, dim): - """Return the solid angle fraction for a given planar angle.""" - alpha = (dim - 1) / 2 - beta = 1/2 - return where(planar_angle < pi/2, - 0.5*betainc(alpha, beta, sin(planar_angle)**2), - 1 - 0.5*betainc(alpha, beta, sin(planar_angle)**2)) - -def solid_angle_fraction2planar_angle(solid_angle_fraction, dim): - """Return the planar angle for a given solid angle fraction.""" - alpha = (dim - 1) / 2 - beta = 1/2 - return where(solid_angle_fraction < 1/2, - arcsin(sqrt(betaincinv(alpha, beta, 2*solid_angle_fraction))), - pi - arcsin(sqrt(betaincinv(alpha, beta, 2*(1-solid_angle_fraction))))) - def rotate_from_nth_canonical(x, axis): """Rotate vector from around the nth canonical axis to the given axis. """ @@ -56,12 +39,6 @@ def rotate_from_nth_canonical(x, axis): - axisn*(xn*s + a*(axisn - 1))/b return x -def random_planar_angle_cdf(maximum_planar_angle, dim): - """Return a random planar angle using inverse transform sampling.""" - return solid_angle_fraction2planar_angle( - 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 @@ -82,19 +59,9 @@ disk.""" x[-1] = cos(planar_angle) return rotate_from_nth_canonical(x, axis) -def random_vector_on_spherical_cap_cdf(axis, maximum_planar_angle): - """Return a random vector uniformly distributed on a spherical -cap. The random planar angle is generated using inverse transform -sampling.""" - return random_on_disk(axis, random_planar_angle_cdf(maximum_planar_angle, axis.size)) - -def random_vector_on_spherical_cap_pdf(axis, maximum_planar_angle): +def random_on_cap(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_on_disk(axis, random_planar_angle_pdf(maximum_planar_angle, axis.size)) -- cgit v1.2.3