Added unifrom sampling option for sphere and torus (#28)

Added uniform sampling option for sphere and torus

Added uniform sampling option for sphere and torus

added uniform sampling option

Added uniform option for sphere and torus

Added uniform option for sphere and Torus

added example usage with uniform option

added uniformity test

ran ruff format

added scikit-learn dependency to pyproject.toml

added example usage with uniform option

added uniformity test

ran ruff format

corrected the test_shapes.py

corrected the dependency

ran ruff format and updated version number

removed the genus g surface test

Author: Sai271828

docs/notebooks/Examples.ipynb

@@ -42,7 +42,7 @@ keywords = [
 ]
 
 [project.optional-dependencies]
-testing = ["pytest", "pytest-cov", "scipy"]
+testing = ["pytest", "pytest-cov", "scikit-learn", "scipy"]
 
 docs = ["sktda_docs_config"]
 

pyproject.toml

@@ -1 +1 @@
-__version__ = "0.2.1"
+__version__ = "0.2.2"

tadasets/_version.py

@@ -66,6 +66,7 @@ def sphere(
     noise: Optional[float] = None,
     ambient: Optional[int] = None,
     seed: Optional[int] = None,
+    uniform: bool = False,
 ) -> np.ndarray:
     """
         Sample ``n`` data points on a sphere.
@@ -82,6 +83,8 @@ def sphere(
         Embed the sphere into a space with ambient dimension equal to `ambient`. The sphere is randomly rotated in this high dimensional space.
     seed : int, optional
         Seed for random state.
+    uniform : bool, default=False
+        If True, sample points uniformly on the sphere. If False, sample points by choosing spherical coordinates uniformly.
 
     Returns
     -------
@@ -90,15 +93,21 @@ def sphere(
     """
 
     rng = np.random.default_rng(seed)
-    theta = rng.random(n) * 2.0 * np.pi
-    phi = rng.random(n) * np.pi
-    rad = np.ones((n,)) * r
 
-    data = np.zeros((n, 3))
+    if uniform:
+        data = rng.standard_normal((n, 3))
+        data = r * data / np.sqrt(np.sum(data**2, 1)[:, None])
+
+    else:
+        theta = rng.random(n) * 2.0 * np.pi
+        phi = rng.random(n) * np.pi
+        rad = np.ones((n,)) * r
+
+        data = np.zeros((n, 3))
 
-    data[:, 0] = rad * np.cos(theta) * np.cos(phi)
-    data[:, 1] = rad * np.cos(theta) * np.sin(phi)
-    data[:, 2] = rad * np.sin(theta)
+        data[:, 0] = rad * np.cos(theta) * np.cos(phi)
+        data[:, 1] = rad * np.cos(theta) * np.sin(phi)
+        data[:, 2] = rad * np.sin(theta)
 
     if noise:
         data += noise * rng.standard_normal(data.shape)
@@ -116,6 +125,7 @@ def torus(
     noise: Optional[float] = None,
     ambient: Optional[int] = None,
     seed: Optional[int] = None,
+    uniform: bool = False,
 ) -> np.ndarray:
     """
     Sample ``n`` data points on a torus.
@@ -134,6 +144,8 @@ def torus(
         Embed the torus into a space with ambient dimension equal to `ambient`. The torus is randomly rotated in this high dimensional space.
     seed : int, optional
         Seed for random state.
+    uniform : bool, default=False
+        If True, sample points uniformly on the torus. If False, sample points by choosing angles uniformly.
 
     Returns
     -------
@@ -144,13 +156,31 @@ def torus(
     assert a <= c, "That's not a torus"
 
     rng = np.random.default_rng(seed)
-    theta = rng.random(n) * 2.0 * np.pi
-    phi = rng.random(n) * 2.0 * np.pi
-
-    data = np.zeros((n, 3))
-    data[:, 0] = (c + a * np.cos(theta)) * np.cos(phi)
-    data[:, 1] = (c + a * np.cos(theta)) * np.sin(phi)
-    data[:, 2] = a * np.sin(theta)
+    if uniform:
+        samples = []
+        while len(samples) < n:
+            u = np.random.uniform(0, 2 * np.pi)
+            v = np.random.uniform(0, 2 * np.pi)
+            x = (c + a * np.cos(v)) * np.cos(u)
+            y = (c + a * np.cos(v)) * np.sin(u)
+            z = a * np.sin(v)
+
+            # REJECTION SAMPLING TO ENSURE UNIFORMITY
+            # The map from u,v to x,y,z is not area-preserving, so we use rejection sampling to ensure uniformity
+            jacobian = a * (c + a * np.cos(v))
+            acceptance_prob = jacobian / (a * (c + a))
+            if np.random.uniform(0, 1) < acceptance_prob:
+                samples.append((x, y, z))
+        data = np.array(samples)
+
+    else:
+        theta = rng.random(n) * 2.0 * np.pi
+        phi = rng.random(n) * 2.0 * np.pi
+
+        data = np.zeros((n, 3))
+        data[:, 0] = (c + a * np.cos(theta)) * np.cos(phi)
+        data[:, 1] = (c + a * np.cos(theta)) * np.sin(phi)
+        data[:, 2] = a * np.sin(theta)
 
     if noise:
         data += noise * rng.standard_normal(data.shape)

tadasets/shapes.py

@@ -2,6 +2,7 @@
 import pytest
 import tadasets
 from scipy.spatial.distance import pdist
+from sklearn.neighbors import NearestNeighbors
 
 
 def norm(p):
@@ -53,6 +54,26 @@ def test_noise(self):
         rs = np.fromiter((norm(p) for p in s), np.float64)
         assert np.any(rs != 1.0)
 
+    def test_uniform(self):
+        s = tadasets.sphere(n=6400, r=1, uniform=True)
+        epsilon = 0.5
+
+        nbrs = NearestNeighbors(radius=epsilon).fit(s)
+
+        # Compute the number of points within the epsilon-ball around each point
+        distances, indices = nbrs.radius_neighbors(s)
+        num_points_in_ball = [len(ind) for ind in indices]
+
+        # Caculate the mean and standard deviation of the counts
+        # For a uniform sampling on the sphere, this ditribution should be roughly normal
+        # The expected mean is n/16 = 400
+        # The standard deviation should be approximately np.sqrt(400) = 20
+        mean_count = np.mean(num_points_in_ball)
+        std_count = np.std(num_points_in_ball)
+
+        assert mean_count > 395 and mean_count < 405
+        assert std_count < 45
+
 
 class TestDsphere:
     def test_d(self):
@@ -100,6 +121,27 @@ def test_ambient(self):
         s = tadasets.torus(n=200, c=3, ambient=15)
         assert s.shape == (200, 15)
 
+    def test_uniform(self):
+        c, a = 3.183, 1
+        s = tadasets.torus(n=16000, c=c, a=a, uniform=True)
+        epsilon = 0.5
+
+        nbrs = NearestNeighbors(radius=epsilon).fit(s)
+
+        # Compute the number of points within the epsilon-ball around each point
+        distances, indices = nbrs.radius_neighbors(s)
+        num_points_in_ball = [len(ind) for ind in indices]
+
+        # Caculate the mean and standard deviation of the counts
+        # For a uniform sampling on the sphere, this ditribution should be roughly normal
+        # The expected mean is n/160 = 100
+        # The standard deviation should be approximately np.sqrt(100) = 10
+        mean_count = np.mean(num_points_in_ball)
+        std_count = np.std(num_points_in_ball)
+
+        assert mean_count > 95 and mean_count < 105
+        assert std_count < 25
+
 
 class TestSwissRoll:
     def test_n(self):