use the dot-product identity to compute angles between vectors (#47)

* use the dot-product identity for angles between vectors

* use the new function to compute the angular distances

* rename a test

* clip to the range of arccos

healpix_convolution/distances.py

@@ -16,6 +16,16 @@ def cell_ids2vectors(cell_ids, nside, nest):
     return np.reshape(vecs, cell_ids.shape + (3,))
 
 
+def angle_between_vectors(a, b, axis):
+    length_a = np.linalg.norm(a, axis=axis)
+    length_b = np.linalg.norm(b, axis=axis)
+    dot_product = np.sum(a * b, axis=axis)
+
+    argument = np.clip(dot_product / (length_a * length_b), -1, 1)
+
+    return np.arccos(argument)
+
+
 def _distances(a, b, axis, nside, nest):
     vec_a = cell_ids2vectors(a, nside, nest)
 
@@ -24,10 +34,7 @@ def _distances(a, b, axis, nside, nest):
     vec_b_ = cell_ids2vectors(np.where(mask, b, 0), nside, nest)
     vec_b = np.where(mask[..., None], vec_b_, np.nan)
 
-    dot_product = np.abs(np.sum(vec_a * vec_b, axis=axis))
-    cross_product = np.linalg.norm(np.cross(vec_a, vec_b, axis=axis), axis=axis)
-
-    return np.arctan2(cross_product, dot_product)
+    return angle_between_vectors(vec_a, vec_b, axis=axis)
 
 
 def angular_distances(neighbours, *, resolution, indexing_scheme="nested", axis=None):

healpix_convolution/tests/test_distances.py

@@ -4,14 +4,24 @@
 import healpix_convolution.distances as hds
 
 
+def test_angle_between_vectors():
+    angles = np.linspace(0, 2 * np.pi, 10)
+    vectors = np.stack([np.cos(angles), np.sin(angles), np.zeros_like(angles)], axis=-1)
+
+    actual = hds.angle_between_vectors(vectors[:1, :], vectors, axis=-1)
+    expected = np.pi - abs(angles - np.pi)
+
+    np.testing.assert_allclose(actual, expected)
+
+
 @pytest.mark.parametrize(
     ["neighbours", "expected"],
     (
         (np.array([[1, 0, 2, 3]]), np.array([[0, 0.25637566, 0.3699723, 0.25574741]])),
         (np.array([[2, 71, 8]]), np.array([[0, 0.25637566, 0.25574741]])),
     ),
 )
-def test_distances_numpy(neighbours, expected):
+def test_angular_distances_numpy(neighbours, expected):
     actual = hds.angular_distances(neighbours, resolution=2, indexing_scheme="nested")
 
     np.testing.assert_allclose(actual, expected)