Add tests.

* tests/test_sambal.py: New file.

1 files changed, 59 insertions, 0 deletions

diff --git a/tests/test_sambal.py b/tests/test_sambal.py
new file mode 100644
index 0000000..e99dd78
--- /dev/null
+++ b/tests/test_sambal.py
@@ -0,0 +1,59 @@
+# sambal --- Sample balls, spheres, spherical caps
+# Copyright © 2021 Arun I <arunisaac@systemreboot.net>
+# Copyright © 2021 Murugesan Venkatapathi <murugesh@iisc.ac.in>
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+# General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see
+# <https://www.gnu.org/licenses/>.
+
+import pytest
+from numpy import arccos, dot, ones, pi, sin, where
+from numpy.linalg import norm
+from numpy.random import default_rng
+from scipy.special import betainc
+from scipy.stats import kstest
+from sambal import random_on_cap
+
+# Set seed of random number generator.
+rng = default_rng(0)
+
+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 make_uniform_cdf(maximum_planar_angle, dim):
+    """Return the CDF of theta uniformly distributed on the spherical
+cap.
+
+    """
+    def cdf(theta):
+        return where(theta > maximum_planar_angle, 1,
+                     planar_angle2solid_angle_fraction(theta, dim)
+                     / planar_angle2solid_angle_fraction(maximum_planar_angle, dim))
+    return cdf
+
+dimensions = [10, 100, 1000, 5000]
+testdata = [*[(dim, 0.35*pi) for dim in dimensions],
+            *[(dim, 0.65*pi) for dim in dimensions]]
+
+@pytest.mark.parametrize("dim,maximum_planar_angle", testdata)
+def test_random_on_cap(dim, maximum_planar_angle):
+    axis = ones(dim)
+    axis = axis / norm(axis)
+    thetas = [arccos(dot(random_on_cap(axis, maximum_planar_angle, rng), axis))
+              for i in range(1000)]
+    assert kstest(thetas, make_uniform_cdf(maximum_planar_angle, dim)).statistic < 0.05