Add phi_units option

tests/configs/ggg_logsas.yaml

@@ -14,8 +14,9 @@ min_sep: 11.
 max_sep: 15.
 nbins: 2
 sep_units: arcmin
-min_phi: 1.
-max_phi: 1.5
+min_phi: 45
+max_phi: 90
+phi_units: deg
 nphi_bins: 10
 bin_type: LogSAS
 

tests/configs/kkk_direct_logsas.yaml

@@ -14,6 +14,7 @@ min_sep: 1.
 max_sep: 10.
 nbins: 10
 nphi_bins: 10
+phi_units: 'radians'
 sep_units: arcmin
 bin_slop: 0
 bin_type: LogSAS

tests/configs/kkk_direct_spherical_logsas.yaml

@@ -15,6 +15,7 @@ max_sep: 100.
 nbins: 3
 sep_units: degrees
 nphi_bins: 6
+phi_units: degrees
 bin_slop: 0
 bin_type: LogSAS
 

tests/configs/kkk_logsas.yaml

@@ -13,8 +13,9 @@ min_sep: 10.
 max_sep: 13.
 nbins: 3
 sep_units: arcmin
-min_phi: 1.0
-max_phi: 1.5
+min_phi: 45.
+max_phi: 90.
+phi_units: degree
 nphi_bins: 5
 bin_type: LogSAS
 

tests/test_ggg.py

@@ -3216,14 +3216,14 @@ def test_ggg_logsas():
     min_sep = 11.
     max_sep = 15.
     nbins = 2
-    min_phi = 1.
-    max_phi = 1.5
+    min_phi = 45
+    max_phi = 90
     nphi_bins = 10
 
     cat = treecorr.Catalog(x=x, y=y, g1=g1, g2=g2, x_units='arcmin', y_units='arcmin')
     ggg = treecorr.GGGCorrelation(min_sep=min_sep, max_sep=max_sep, nbins=nbins,
                                   min_phi=min_phi, max_phi=max_phi, nphi_bins=nphi_bins,
-                                  sep_units='arcmin', bin_type='LogSAS')
+                                  sep_units='arcmin', phi_units='degrees', bin_type='LogSAS')
     t0 = time.time()
     ggg.process(cat)
     t1 = time.time()
@@ -3233,13 +3233,11 @@ def test_ggg_logsas():
     # basically the same answer.
     gggc = treecorr.GGGCorrelation(min_sep=min_sep, max_sep=max_sep, nbins=nbins,
                                    min_phi=min_phi, max_phi=max_phi, nphi_bins=nphi_bins,
-                                   sep_units='arcmin', bin_type='LogSAS')
+                                   sep_units='arcmin', phi_units='degrees', bin_type='LogSAS')
     t0 = time.time()
     gggc.process(cat,cat,cat, ordered=True)
     t1 = time.time()
     print('cross process time = ',t1-t0)
-    print(ggg.gam0)
-    print(gggc.gam0)
     np.testing.assert_allclose(gggc.ntri, ggg.ntri, rtol=1.e-3)
     np.testing.assert_allclose(gggc.meanlogd1, ggg.meanlogd1, rtol=1.e-3)
     np.testing.assert_allclose(gggc.meanlogd2, ggg.meanlogd2, rtol=1.e-3)
@@ -3279,7 +3277,7 @@ def test_ggg_logsas():
     s = d2
     # And let t be from c1 to c2. t = |t| e^Iphi
     # |t| = d3
-    t = d3 * np.exp(1j * ggg.meanphi)
+    t = d3 * np.exp(1j * ggg.meanphi * np.pi/180.)
 
     q1 = (s + t)/3.
     q2 = q1 - t
@@ -3455,7 +3453,7 @@ def test_ggg_logsas():
         # The read function should reshape them to the right shape.
         ggg2 = treecorr.GGGCorrelation(min_sep=min_sep, max_sep=max_sep, nbins=nbins,
                                        min_phi=min_phi, max_phi=max_phi, nphi_bins=nphi_bins,
-                                       sep_units='arcmin', bin_type='LogSAS')
+                                       sep_units='arcmin', phi_units='degrees', bin_type='LogSAS')
         ggg2.read(out_file_name1)
         np.testing.assert_almost_equal(ggg2.logd2, ggg.logd2)
         np.testing.assert_almost_equal(ggg2.logd3, ggg.logd3)

tests/test_kkk.py

@@ -1201,7 +1201,8 @@ def test_direct_logsas():
     nbins = 10
     nphi_bins = 10
     kkk = treecorr.KKKCorrelation(min_sep=min_sep, max_sep=max_sep, nbins=nbins,
-                                  nphi_bins=nphi_bins, brute=True, bin_type='LogSAS')
+                                  nphi_bins=nphi_bins, phi_units='rad',
+                                  brute=True, bin_type='LogSAS')
     kkk.process(cat, num_threads=2)
 
     log_min_sep = np.log(min_sep)
@@ -1374,7 +1375,8 @@ def test_direct_logsas_spherical():
     nbins = 3
     nphi_bins = 6
     kkk = treecorr.KKKCorrelation(min_sep=min_sep, max_sep=max_sep, sep_units='deg',
-                                  nbins=nbins, nphi_bins=nphi_bins, brute=True, bin_type='LogSAS')
+                                  nbins=nbins, nphi_bins=nphi_bins, phi_units='deg',
+                                  brute=True, bin_type='LogSAS')
     kkk.process(cat)
 
     r = np.sqrt(x**2 + y**2 + z**2)
@@ -1482,16 +1484,16 @@ def test_direct_logsas_spherical():
     true_meand3_arc[posa] /= true_weight_arc[posa]
     true_meanphi_arc[posa] /= true_weight_arc[posa]
 
-    # Convert chord distances and angle to spherical values (in degrees for distances).
+    # Convert chord distances and angle to spherical values (in degrees)
     # cosphi = (d2^2 + d3^2 - d1^2 - 1/2 d2^2 d3^2) / (2 d2 d3 sqrt(1-d2^2) sqrt(1-d3^2))
     # Fix this first, while the ds are still chord distances.
     cosphi = np.cos(true_meanphi)
     cosphi -= 0.25 * true_meand2 * true_meand3
     cosphi /= np.sqrt(1-0.25*true_meand2**2) * np.sqrt(1-0.25*true_meand3**2)
-    true_meanphi[:] = np.arccos(cosphi)
+    true_meanphi[:] = np.arccos(cosphi) * 180./np.pi
     for dd in [true_meand1, true_meand2, true_meand3]:
         dd[:] = 2 * np.arcsin(dd/2) * 180/np.pi
-    for dd in [true_meand1_arc, true_meand2_arc, true_meand3_arc]:
+    for dd in [true_meand1_arc, true_meand2_arc, true_meand3_arc, true_meanphi_arc]:
         dd *= 180. / np.pi
 
     np.testing.assert_array_equal(kkk.ntri, true_ntri)
@@ -1522,7 +1524,7 @@ def test_direct_logsas_spherical():
     # Repeat with binslop = 0
     # And don't do any top-level recursion so we actually test not going to the leaves.
     kkk = treecorr.KKKCorrelation(min_sep=min_sep, max_sep=max_sep, sep_units='deg',
-                                  nbins=nbins, nphi_bins=nphi_bins,
+                                  nbins=nbins, nphi_bins=nphi_bins, phi_units='deg',
                                   bin_slop=0, max_top=0, bin_type='LogSAS')
     kkk.process(cat)
     np.testing.assert_array_equal(kkk.ntri, true_ntri)
@@ -1535,7 +1537,7 @@ def test_direct_logsas_spherical():
 
     # Now do Arc metric, where distances and angles use spherical geometry.
     kkk = treecorr.KKKCorrelation(min_sep=min_sep, max_sep=max_sep, sep_units='deg',
-                                  nbins=nbins, nphi_bins=nphi_bins,
+                                  nbins=nbins, nphi_bins=nphi_bins, phi_units='deg',
                                   bin_slop=0, bin_type='LogSAS', metric='Arc')
     kkk.process(cat, num_threads=1)
     np.testing.assert_array_equal(kkk.ntri, true_ntri_arc)
@@ -1547,7 +1549,7 @@ def test_direct_logsas_spherical():
     np.testing.assert_allclose(kkk.meanphi[posa], true_meanphi_arc[posa], rtol=1.e-4, atol=1.e-6)
 
     kkk = treecorr.KKKCorrelation(min_sep=min_sep, max_sep=max_sep, sep_units='deg',
-                                  nbins=nbins, nphi_bins=nphi_bins,
+                                  nbins=nbins, nphi_bins=nphi_bins, phi_units='deg',
                                   bin_slop=0, max_top=0, bin_type='LogSAS', metric='Arc')
     kkk.process(cat)
     np.testing.assert_array_equal(kkk.ntri, true_ntri_arc)
@@ -1974,14 +1976,14 @@ def test_kkk_logsas():
     min_sep = 10.
     max_sep = 13.
     nbins = 3
-    min_phi = 1.
-    max_phi = 1.5
+    min_phi = 45.
+    max_phi = 90.
     nphi_bins = 5
 
     cat = treecorr.Catalog(x=x, y=y, k=kappa, x_units='arcmin', y_units='arcmin')
     kkk = treecorr.KKKCorrelation(min_sep=min_sep, max_sep=max_sep, nbins=nbins,
                                   min_phi=min_phi, max_phi=max_phi, nphi_bins=nphi_bins,
-                                  sep_units='arcmin', bin_type='LogSAS')
+                                  sep_units='arcmin', phi_units='deg', bin_type='LogSAS')
     t0 = time.time()
     kkk.process(cat)
     print(kkk.ntri)
@@ -1992,7 +1994,7 @@ def test_kkk_logsas():
     # basically the same answer.
     kkkc = treecorr.KKKCorrelation(min_sep=min_sep, max_sep=max_sep, nbins=nbins,
                                    min_phi=min_phi, max_phi=max_phi, nphi_bins=nphi_bins,
-                                   sep_units='arcmin', bin_type='LogSAS')
+                                   sep_units='arcmin', phi_units='deg', bin_type='LogSAS')
     t0 = time.time()
     kkkc.process(cat,cat,cat, ordered=True)
     t1 = time.time()
@@ -2072,7 +2074,7 @@ def test_kkk_logsas():
         # The read function should reshape them to the right shape.
         kkk2 = treecorr.KKKCorrelation(min_sep=min_sep, max_sep=max_sep, nbins=nbins,
                                        min_phi=min_phi, max_phi=max_phi, nphi_bins=nphi_bins,
-                                       sep_units='arcmin', bin_type='LogSAS')
+                                       sep_units='arcmin', phi_units='deg', bin_type='LogSAS')
         kkk2.read(out_file_name)
         np.testing.assert_almost_equal(kkk2.logd2, kkk.logd2)
         np.testing.assert_almost_equal(kkk2.logd3, kkk.logd3)
@@ -2091,6 +2093,7 @@ def test_kkk_logsas():
         assert kkk2.coords == kkk.coords
         assert kkk2.metric == kkk.metric
         assert kkk2.sep_units == kkk.sep_units
+        assert kkk2.phi_units == kkk.phi_units
         assert kkk2.bin_type == kkk.bin_type
 
 if __name__ == '__main__':

treecorr/corr3base.py

@@ -157,6 +157,10 @@ class Corr3(object):
         phi_bin_size (float): Analogous to bin_size for the phi values. (default: bin_size)
         min_phi (float):    Analogous to min_sep for the phi values. (default: 0)
         max_phi (float):    Analogous to max_sep for the phi values. (default: np.pi)
+        phi_units (str):    The units to use for the phi values, given as a string.  This
+                            includes both min_phi and max_phi above, as well as the units of the
+                            output meanphi values.  Valid options are arcsec, arcmin, degrees,
+                            hours, radians.  (default: radians)
 
         nubins (int):       Analogous to nbins for the u values when bin_type=LogRUV.  (The default
                             is to calculate from ubin_size = bin_size, min_u = 0, max_u = 1, but
@@ -288,6 +292,8 @@ class Corr3(object):
                 'The minimum phi to include in the output.'),
         'max_phi' : (float, False, None, None,
                 'The maximum phi to include in the output.'),
+        'phi_units' : (str, False, None, coord.AngleUnit.valid_names,
+                'The units to use for min_phi and max_phi.'),
         'brute' : (bool, False, False, [False, True],
                 'Whether to use brute-force algorithm'),
         'verbose' : (int, False, 1, [0, 1, 2, 3],
@@ -472,10 +478,15 @@ def __init__(self, config=None, *, logger=None, rng=None, **kwargs):
                     raise TypeError("%s is invalid for bin_type=LogSAS"%key)
             self._ro._bintype = _treecorr.LogSAS
             # Note: we refer to phi as u in the _ro namespace to make function calls easier.
+            self._ro.phi_units = self.config.get('phi_units','')
+            self._ro._phi_units = get(self.config,'phi_units',str,'radians')
             self._ro.min_u = float(self.config.get('min_phi', 0.))
-            self._ro.max_u = float(self.config.get('max_phi', np.pi))
+            self._ro.max_u = float(self.config.get('max_phi', np.pi / self._phi_units))
+            self._ro.min_u = self.min_phi * self._phi_units
+            self._ro.max_u = self.max_phi * self._phi_units
             if self.min_phi < 0 or self.max_phi > np.pi:
-                raise ValueError("Invalid range for phi: %f - %f"%(self.min_phi, self.max_phi))
+                raise ValueError("Invalid range for phi: %f - %f"%(
+                                 self.min_phi/self._phi_units, self.max_phi/self._phi_units))
             if self.min_phi >= self.max_phi:
                 raise ValueError("max_phi must be larger than min_phi")
             self._ro.ubin_size = float(self.config.get('phi_bin_size', bin_size))
@@ -761,6 +772,10 @@ def phi1d(self):
     def phi(self): 
         assert self.bin_type == 'LogSAS'
         return self._ro.phi
+    @property
+    def phi_units(self): return self._ro.phi_units
+    @property
+    def _phi_units(self): return self._ro._phi_units
 
     def _equal_binning(self, other, brief=False):
         # A helper function to test if two Corr3 objects have the same binning parameters
@@ -788,6 +803,7 @@ def _equal_binning(self, other, brief=False):
             return eq
         else:
             return (self.sep_units == other.sep_units and
+                    (self.bin_type == 'LogRUV' or self.phi_units == other.phi_units) and
                     self.coords == other.coords and
                     self.bin_slop == other.bin_slop and
                     self.xperiod == other.xperiod and
@@ -1429,6 +1445,8 @@ def _apply_units(self, mask):
             self.meand3[mask] = 2. * np.arcsin(self.meand3[mask]/2.)
             self.meanlogd3[mask] = np.log(2.*np.arcsin(np.exp(self.meanlogd3[mask])/2.))
 
+        if self.bin_type == 'LogSAS':
+            self.meanphi[mask] /= self._phi_units
         self.meand1[mask] /= self._sep_units
         self.meanlogd1[mask] -= self._log_sep_units
         self.meand2[mask] /= self._sep_units

treecorr/gggcorrelation.py

@@ -1024,7 +1024,7 @@ def calculateMap3(self, *, R=None, k2=1, k3=1):
             t = tx + 1j * ty
         else:
             d3 = np.outer(1./R, self.meand3.ravel())
-            t = d3 * np.exp(1j * self.meanphi.ravel())
+            t = d3 * np.exp(1j * self.meanphi.ravel() * self._phi_units)
 
         # Next we need to construct the T values.
         T0, T1, T2, T3 = self._calculateT(s,t,1.,k2,k3)