aboutsummaryrefslogtreecommitdiff
path: root/contrib
diff options
context:
space:
mode:
authorArun Isaac2021-03-19 21:50:55 +0530
committerArun Isaac2021-03-19 21:57:31 +0530
commite4ca0edf037b924b217bd5c2223e570a8e888ee9 (patch)
tree513efa33557bb18da5b82927e9348683a4c18888 /contrib
parent5b8aa0c2cb6f187eba957de9400913b83bbf539a (diff)
downloadnsmc-e4ca0edf037b924b217bd5c2223e570a8e888ee9.tar.gz
nsmc-e4ca0edf037b924b217bd5c2223e570a8e888ee9.tar.lz
nsmc-e4ca0edf037b924b217bd5c2223e570a8e888ee9.zip
Implement rejection sampling based cone sampling.
* contrib/cone-vector.py: Import log. (random_planar_angle_pdf, random_vector_on_spherical_cap_pdf): New functions.
Diffstat (limited to 'contrib')
-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