Merge pull request #35 from chrisjonesBSU/bond_dist_range

Add radians option to `angle_between_vectors` and add value ranges to bond distribution functions.

cmeutils/geometry.py

@@ -60,8 +60,8 @@ def get_plane_normal(points):
     return ctr, normal
 
 
-def angle_between_vectors(u, v, min_angle=True):
-    """Calculate the angle between two vectors in degrees.
+def angle_between_vectors(u, v, min_angle=True, degrees=True):
+    """Calculate the angle between two vectors in radians or degrees.
 
     Parameters
     ----------
@@ -73,15 +73,21 @@ def angle_between_vectors(u, v, min_angle=True):
         Whether to return the supplement if the angle is greater than 90
         degrees. Useful for calculating the minimum angle between the normal
         vectors of planes as direction doesn't matter.
+    degrees : bool, default True
+        If True, the angle is returned in degrees.
+        If False, the angle is returned in radians.
 
     Returns
     -------
     angle: float
-        Angle between u and v in degrees
+        Angle between u and v
+
     """
-    angle = np.rad2deg(
-        np.arccos(u.dot(v) / (np.linalg.norm(u) * np.linalg.norm(v)))
-    )
-    if angle > 90 and min_angle:
-        return 180 - angle
+    # Angle in radians
+    angle = np.arccos(u.dot(v) / (np.linalg.norm(u) * np.linalg.norm(v)))
+    if angle > np.pi/2 and min_angle:
+        angle = np.pi - angle
+
+    if degrees:
+        return np.rad2deg(angle)
     return angle

cmeutils/plotting.py

@@ -1,7 +1,8 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
-def get_histogram(data, normalize=False, bins="auto"):
+
+def get_histogram(data, normalize=False, bins="auto", x_range=None):
     """Bins a 1-D array of data into a histogram using
     the numpy.histogram method.
 
@@ -16,6 +17,9 @@ def get_histogram(data, normalize=False, bins="auto"):
     bins : float, int, or str, default="auto"
         Method used by numpy to determine bin borders.
         Check the numpy.histogram docs for more details.
+    x_range : (float, float), default = None
+        The lower and upper range of the histogram bins.
+        If set to None, then the min and max values of data are used.
 
     Returns
     -------
@@ -25,7 +29,7 @@ def get_histogram(data, normalize=False, bins="auto"):
         Array of the bin height values
 
     """
-    bin_heights, bin_borders = np.histogram(data, bins=bins)
+    bin_heights, bin_borders = np.histogram(data, bins=bins, range=x_range)
     if normalize is True:
         bin_heights = bin_heights/sum(bin_heights)
     bin_widths = np.diff(bin_borders)

cmeutils/structure.py

@@ -1,3 +1,5 @@
+import warnings
+
 import freud
 import gsd
 import gsd.hoomd
@@ -16,7 +18,10 @@ def angle_distribution(
         C_name,
         start=0,
         stop=-1,
+        degrees=False,
         histogram=False,
+        theta_min=0.0,
+        theta_max=None,
         normalize=False,
         bins="auto"
 ):
@@ -35,10 +40,19 @@ def angle_distribution(
         Negative numbers index from the end. (default 0)
     stop : int
         Final frame index for accumulating bond lengths. (default -1)
+    degrees : bool, default=False
+        If True, the angle values are returned in degrees.
+        if False, the angle values are returned in radians.
     histogram : bool, default=False
         If set to True, places the resulting angles into a histogram
         and retrums the histogram's bin centers and heights as 
         opposed to the actual calcualted angles.
+    theta_min : float, default = 0.0
+        Sets the minimum theta value to be included in the distribution
+    theta_max : float, default = None 
+        Sets the maximum theta value to be included in the distribution
+        If left as None, then theta_max will be either pi radians or
+        180 degrees depending on the value set for the degrees parameter
     normalize : bool, default=False
         If set to True, normalizes the angle distribution by the
         sum of the bin heights, so that the distribution adds up to 1. 
@@ -55,7 +69,11 @@ def angle_distribution(
         If histogram is True, returns a 2D array of bin centers and bin heights.
 
     """
-    angles = []
+    if not degrees and theta_max is None:
+        theta_max = np.pi
+    elif degrees and theta_max is None:
+        theta_max = 180
+
     trajectory = gsd.hoomd.open(gsd_file, mode="rb")
     name = "-".join([A_name, B_name, C_name])
     name_rev = "-".join([C_name, B_name, A_name])
@@ -81,13 +99,23 @@ def angle_distribution(
                 pos1_unwrap = pos1 + (img1 * snap.configuration.box[:3])
                 pos2_unwrap = pos2 + (img2 * snap.configuration.box[:3])
                 pos3_unwrap = pos3 + (img3 * snap.configuration.box[:3])
-                u = pos2_unwrap - pos1_unwrap
+                u = pos1_unwrap - pos2_unwrap
                 v = pos3_unwrap - pos2_unwrap
-                angles.append(np.round(angle_between_vectors(u, v, False), 3))
+                angles.append(
+                        np.round(angle_between_vectors(u, v, False, degrees), 3)
+                )
     trajectory.close()
 
     if histogram:
-        bin_centers, bin_heights = get_histogram(np.array(angles), bins=bins)
+        if angles[0] < theta_min or angles[-1] > theta_max:
+            warnings.warn("There are bond angles that fall outside of "
+                    "your set theta_min and theta_max range. "
+                    "You may want to adjust this range to "
+                    "include all bond angles."
+            )
+        bin_centers, bin_heights = get_histogram(
+                np.array(angles), bins=bins, x_range=(theta_min, theta_max)
+        )
         return np.stack((bin_centers, bin_heights)).T
     else:
         return np.array(angles)
@@ -100,6 +128,8 @@ def bond_distribution(
     start=0,
     stop=-1,
     histogram=False,
+    l_min=0.0,
+    l_max=4.0,
     normalize=True,
     bins=100
 ):
@@ -121,6 +151,10 @@ def bond_distribution(
         If set to True, places the resulting bonds into a histogram
         and retrums the histogram's bin centers and heights as 
         opposed to the actual calcualted bonds.
+    l_min : float, default = 0.0
+        Sets the minimum bond length to be included in the distribution
+    l_max : float, default = 5.0 
+        Sets the maximum bond length value to be included in the distribution
     normalize : bool, default=False
         If set to True, normalizes the angle distribution by the
         sum of the bin heights, so that the distribution adds up to 1. 
@@ -162,7 +196,14 @@ def bond_distribution(
     trajectory.close()
 
     if histogram:
-        bin_centers, bin_heights = get_histogram(np.array(bonds), bins=bins)
+        if bonds[0] < l_min or bonds[-1] > l_max:
+            warnings.warn("There are bond lengths that fall outside of "
+                    "your set l_min and l_max range. You may want to adjust "
+                    "this range to include all bond lengths."
+            )
+        bin_centers, bin_heights = get_histogram(
+                np.array(bonds), bins=bins, x_range=(l_min, l_max)
+        )
         return np.stack((bin_centers, bin_heights)).T
     else:
         return np.array(bonds)

cmeutils/tests/test_geometry.py

@@ -1,3 +1,5 @@
+import math
+
 import numpy as np
 import pytest
 
@@ -20,10 +22,25 @@ def test_get_plane_normal(self):
         with pytest.raises(AssertionError):
             get_plane_normal(points[:2])
 
-    def test_angle_between_vectors(self):
+    def test_angle_between_vectors_deg(self):
         assert np.isclose(
             90, angle_between_vectors(np.array([1,0,0]),np.array([0,1,0]))
         )
         assert np.isclose(
             0, angle_between_vectors(np.array([1,0,0]),np.array([-1,0,0]))
         )
+
+    def test_angle_between_vectors_rad(self):
+        assert np.isclose(
+            math.pi/2,
+            angle_between_vectors(
+                np.array([1,0,0]),np.array([0,1,0]), degrees=False
+            )
+        )
+        assert np.isclose(
+            0,
+            angle_between_vectors(
+                np.array([1,0,0]),np.array([-1,0,0]), degrees=False
+            )
+        )
+

cmeutils/tests/test_structure.py

@@ -18,11 +18,32 @@
     )
 
 class TestStructure(BaseTest):
-    def test_angle_distribution(self, p3ht_gsd):
-        angles = angle_distribution(p3ht_gsd, "cc", "ss", "cc", start=0, stop=1)
+    def test_angle_distribution_deg(self, p3ht_gsd):
+        angles = angle_distribution(p3ht_gsd, "cc", "ss", "cc", start=0, stop=1, degrees=True)
         for ang in angles:
             assert 80 < ang < 100
 
+    def test_angle_distribution_range(self, p3ht_gsd):
+        angles = angle_distribution(
+                p3ht_gsd,
+                "cc",
+                "ss",
+                "cc",
+                start=0,
+                stop=1,
+                histogram=True,
+                degrees=True,
+                theta_min=10,
+                theta_max=180
+        )
+        assert np.allclose(angles[0,0],10, atol=0.5)
+        assert np.allclose(angles[-1,0], 180, atol=0.5)
+
+    def test_angle_distribution_rad(self, p3ht_gsd):
+        angles = angle_distribution(p3ht_gsd, "cc", "ss", "cc", start=0, stop=1, degrees=False)
+        for ang in angles:
+            assert 1.40 < ang < 1.75
+
     def test_angle_distribution_order(self, p3ht_gsd):
         angles = angle_distribution(p3ht_gsd, "ss", "cc", "cd", start=0, stop=1)
         angles2 = angle_distribution(p3ht_gsd, "cd", "cc", "ss", start=0, stop=1)
@@ -37,8 +58,23 @@ def test_angle_not_found(self, p3ht_gsd):
     def test_bond_distribution(self, p3ht_gsd):
         bonds = bond_distribution(p3ht_gsd, "cc", "ss", start=0, stop=1)
         for bond in bonds:
+            print(bond)
             assert 0.45 < bond < 0.52
 
+    def test_bond_distribution_range(self, p3ht_gsd):
+        bonds = bond_distribution(
+                p3ht_gsd,
+                "cc",
+                "ss",
+                start=0,
+                stop=1,
+                l_min=0,
+                l_max=1,
+                histogram=True
+        )
+        assert np.allclose(bonds[0,0],0, atol=0.5)
+        assert np.allclose(bonds[-1,0], 1, atol=0.5)
+
     def test_bond_distribution_order(self, p3ht_gsd):
         bonds = bond_distribution(p3ht_gsd, "cc", "ss", start=0, stop=1)
         bonds2 = bond_distribution(p3ht_gsd, "ss", "cc", start=0, stop=1)
@@ -59,6 +95,19 @@ def test_bond_histogram(self, p3ht_gsd):
         assert bonds_hist.ndim == 2
         assert bonds_no_hist.ndim == 1
 
+    def test_bond_range_outside(self, p3ht_gsd):
+        with pytest.warns(UserWarning):
+            bonds_hist = bond_distribution(
+                p3ht_gsd,
+                "cc",
+                "ss",
+                start=0,
+                stop=1,
+                histogram=True,
+                l_min=0.52,
+                l_max=0.60
+            )
+
     def test_angle_histogram(self, p3ht_gsd):
         angles_hist = angle_distribution(
                 p3ht_gsd, "cc", "ss", "cc", start=0, stop=1, histogram=True
@@ -68,6 +117,20 @@ def test_angle_histogram(self, p3ht_gsd):
         )
         assert angles_hist.ndim == 2
         assert angles_no_hist.ndim == 1
+
+    def test_angle_range_outside(self, p3ht_gsd):
+        with pytest.warns(UserWarning):
+            angles_hist = angle_distribution(
+                    p3ht_gsd,
+                    "cc",
+                    "ss",
+                    "cc",
+                    start=0,
+                    stop=1,
+                    histogram=True,
+                    theta_min = 120,
+                    theta_max=180
+            )
 
     def test_gsd_rdf(self, gsdfile_bond):
         rdf_ex, norm = gsd_rdf(gsdfile_bond, "A", "B")