Merge dcd7910 into d6fda00

AegeanTools/angle_tools.py

@@ -145,7 +145,7 @@ def gcd(ra1, dec1, ra2, dec2):
     dlat = dec2 - dec1
     a = np.sin(np.radians(dlat) / 2) ** 2
     a += np.cos(np.radians(dec1)) * np.cos(np.radians(dec2)) * np.sin(np.radians(dlon) / 2) ** 2
-    sep = np.degrees(2 * np.arcsin(min(1, np.sqrt(a))))
+    sep = np.degrees(2 * np.arcsin(np.minimum(1, np.sqrt(a))))
     return sep
 
 

AegeanTools/cluster.py

@@ -45,7 +45,7 @@ def norm_dist(src1, src2):
         The normalised distance.
 
     """
-    if src1 == src2:
+    if np.all(src1 == src2):
         return 0
     dist = gcd(src1.ra, src1.dec, src2.ra, src2.dec) # degrees
 
@@ -81,7 +81,7 @@ def sky_dist(src1, src2):
     :func:`AegeanTools.angle_tools.gcd`
     """
 
-    if src1 == src2:
+    if np.all(src1 == src2):
         return 0
     return gcd(src1.ra, src1.dec, src2.ra, src2.dec) # degrees
 

tests/test_cluster.py

@@ -46,6 +46,34 @@ def test_sky_dist():
         raise AssertionError()
 
 
+def test_vectorized():
+    # random data as struct array with interface like SimpleSource
+    X = np.random.RandomState(0).rand(20, 6)
+    Xr = np.rec.array(X.view([('ra', 'f8'), ('dec', 'f8'),
+                              ('a', 'f8'), ('b', 'f8'),
+                              ('pa', 'f8'),
+                              ('peak_flux', 'f8')]).ravel())
+
+    def to_ss(x):
+        "Convert numpy.rec to SimpleSource"
+        out = SimpleSource()
+        for f in Xr.dtype.names:
+            setattr(out, f, getattr(x, f))
+        return out
+
+    for dist in [cluster.norm_dist, cluster.sky_dist]:
+        x0 = Xr[0]
+        # calculate distance of x0 to all of Xr with vectorized operations:
+        dx0all = dist(x0, Xr)
+        for i in range(len(Xr)):
+            xi = Xr[i]
+            dx0xi = dist(x0, xi)
+            # check equivalence between pairs of sources and vectorized
+            assert np.isclose(dx0xi, dx0all[i], atol=0)
+            # check equivalence between SimpleSource and numpy.record
+            assert np.isclose(dx0xi, dist(to_ss(x0), to_ss(xi)), atol=0)
+
+
 def test_pairwise_elliptical_binary():
     src1 = SimpleSource()
     src1.ra = 0