Merge pull request #181 from TheoChem-VU/180-add-sampling-of-minimum-…

…volume-spheroids-of-molecules

Sampling of minimum volume spheroids of molecules

src/tcutility/geometry.py

@@ -391,7 +391,7 @@ def RMSD(X: np.ndarray, Y: np.ndarray, axis: Union[int, None] = None, use_kabsch
 
 def random_points_on_sphere(shape: Tuple[int], radius: float = 1) -> np.ndarray:
     """
-    Generate a random points on a sphere with a specified radius.
+    Generate random points on a sphere with a specified radius.
 
     Args:
         shape: The shape of the resulting points, generally shape[0] coordinates with shape[1] dimensions
@@ -401,13 +401,14 @@ def random_points_on_sphere(shape: Tuple[int], radius: float = 1) -> np.ndarray:
         Array of coordinates on a sphere.
     """
     x = np.random.randn(*shape)
-    x = x / np.linalg.norm(x, axis=1) * radius
+    x = x / np.linalg.norm(x, axis=1, keepdims=True) * radius
     return x
 
 
 def random_points_in_anular_sphere(shape: Tuple[int], min_radius: float = 0, max_radius: float = 1):
     """
-    Generate a random points in a sphere with specified radii.
+    Generate random points in an sphere or anular sphere with specified radii. 
+    An anular sphere is a hollow sphere of a certain thickness.
 
     Args:
         shape: The shape of the resulting points, generally shape[0] coordinates with shape[1] dimensions
@@ -419,3 +420,35 @@ def random_points_in_anular_sphere(shape: Tuple[int], min_radius: float = 0, max
     """
     random_radii = np.random.rand(shape[0]) * (max_radius - min_radius) + min_radius
     return random_points_on_sphere(shape, random_radii)
+
+
+def random_points_on_spheroid(coordinates: np.ndarray, Nsamples: int = 1, margin: float = 0):
+    """
+    Generate random points on a spheroid generated by a set of coordinates.
+
+    Args:
+        coordinates: The (n x dim) set of coordinates that is used to generate the minimum-volume spheroid.
+        Nsamples: The number of samples to return.
+        margin: the spacing between the sampling spheroid and the minimum-volume spheroid.
+
+    Returns:
+        Array of coordinates on a spheroid.
+    """
+    # for this to work we should first get the centroid of our molecule
+    centroid = np.mean(coordinates, axis=0)
+    # and get the centered coordiantes
+    Xc = coordinates - centroid
+
+    # we then do a singular-value decomposition to obtain
+    # the three principle components (Vh) with their eigenvalues (s)
+    _, s, Vh = scipy.linalg.svd(Xc)
+
+    # then compute a transformation matrix for generating the correct spheroid
+    transform = Transform()
+    transform.translate(centroid)
+    transform.rotate((np.diag(s/2 + margin) @ Vh).T)
+
+    # to sample the spheroid we generate points on a 
+    # sphere and transform them to our spheroid
+    p = random_points_on_sphere((Nsamples, Xc.shape[1]))
+    return transform(p)