From 7e04a30c82cd7242b2de8cff69e982ab9c5688b8 Mon Sep 17 00:00:00 2001 From: Mike Jarvis Date: Thu, 4 Jan 2024 03:05:55 -0500 Subject: [PATCH] Add phi_units option --- tests/configs/ggg_logsas.yaml | 5 ++-- tests/configs/kkk_direct_logsas.yaml | 1 + .../configs/kkk_direct_spherical_logsas.yaml | 1 + tests/configs/kkk_logsas.yaml | 5 ++-- tests/test_ggg.py | 14 ++++----- tests/test_kkk.py | 29 ++++++++++--------- treecorr/corr3base.py | 22 ++++++++++++-- treecorr/gggcorrelation.py | 2 +- 8 files changed, 51 insertions(+), 28 deletions(-) diff --git a/tests/configs/ggg_logsas.yaml b/tests/configs/ggg_logsas.yaml index 93235f95..7572ecb0 100644 --- a/tests/configs/ggg_logsas.yaml +++ b/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 diff --git a/tests/configs/kkk_direct_logsas.yaml b/tests/configs/kkk_direct_logsas.yaml index 4e2dae25..e9adf3da 100644 --- a/tests/configs/kkk_direct_logsas.yaml +++ b/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 diff --git a/tests/configs/kkk_direct_spherical_logsas.yaml b/tests/configs/kkk_direct_spherical_logsas.yaml index e2a8b34a..d0f72bed 100644 --- a/tests/configs/kkk_direct_spherical_logsas.yaml +++ b/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 diff --git a/tests/configs/kkk_logsas.yaml b/tests/configs/kkk_logsas.yaml index d6fb31ca..e00b6a70 100644 --- a/tests/configs/kkk_logsas.yaml +++ b/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 diff --git a/tests/test_ggg.py b/tests/test_ggg.py index 2f84a47a..ebd83afe 100644 --- a/tests/test_ggg.py +++ b/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) diff --git a/tests/test_kkk.py b/tests/test_kkk.py index a4767711..3dccd3ea 100644 --- a/tests/test_kkk.py +++ b/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__': diff --git a/treecorr/corr3base.py b/treecorr/corr3base.py index f7f4067c..8529aca4 100644 --- a/treecorr/corr3base.py +++ b/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 diff --git a/treecorr/gggcorrelation.py b/treecorr/gggcorrelation.py index d15896e6..cfbd11d0 100644 --- a/treecorr/gggcorrelation.py +++ b/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)