From 5ed7d731a7d8b4aeb74fff297734b7111c9c3a4a Mon Sep 17 00:00:00 2001 From: Pierre-Francois Leget Date: Wed, 8 Jan 2025 13:30:00 -0500 Subject: [PATCH 1/5] comp E/B mode using TreeCorr --- requirements.txt | 2 +- treegp/__init__.py | 3 ++- treegp/_version.py | 2 +- treegp/utils.py | 54 ++++++++++++++++++++++++++++++++++++++++++---- 4 files changed, 54 insertions(+), 7 deletions(-) diff --git a/requirements.txt b/requirements.txt index 49cee02..35218ea 100644 --- a/requirements.txt +++ b/requirements.txt @@ -4,6 +4,6 @@ scipy matplotlib iminuit>1.3,<1.4; python_version < '3.0' # iminuit 1.4 fails on python 2.7. iminuit>2; python_version >= '3.0' # iminuit changed API in v2.0. We use the 2.x API. -treecorr>=4.2 +treecorr>=5.0 fitsio>=0.9.12 scikit-learn>=0.18 diff --git a/treegp/__init__.py b/treegp/__init__.py index 5ac47d3..189d2cf 100644 --- a/treegp/__init__.py +++ b/treegp/__init__.py @@ -18,7 +18,7 @@ from .meanify import meanify -from .utils import comp_eb +from .utils import comp_eb, comp_eb_treecorr __all__ = [ "__version__", @@ -32,4 +32,5 @@ "eval_kernel", "meanify", "comp_eb", + "comp_eb_treecorr", ] diff --git a/treegp/_version.py b/treegp/_version.py index 39424a3..1e5422b 100644 --- a/treegp/_version.py +++ b/treegp/_version.py @@ -1,2 +1,2 @@ -__version__ = "1.2.0" +__version__ = "1.3.0" __version_info__ = tuple(map(int, __version__.split("."))) diff --git a/treegp/utils.py b/treegp/utils.py index 04e9826..ead06db 100644 --- a/treegp/utils.py +++ b/treegp/utils.py @@ -1,7 +1,5 @@ import numpy as np - -# TO DO: E/B mode computation using numpy should use -# treecorr instead. +import treecorr def vcorr(x, y, dx, dy, rmin=5.0 / 3600.0, rmax=1.5, dlogr=0.05, maxpts=30000): @@ -37,7 +35,7 @@ def vcorr(x, y, dx, dy, rmin=5.0 / 3600.0, rmax=1.5, dlogr=0.05, maxpts=30000): y = y[use] dx = dx[use] dy = dy[use] - print("Length ", len(x)) + # print("Length ", len(x)) # Get index arrays that make all unique pairs i1, i2 = np.triu_indices(len(x)) # Omit self-pairs @@ -109,3 +107,51 @@ def comp_eb(u, v, du, dv, **kwargs): xib = xiB(logr, xiplus, ximinus) xie = xiplus - xib return xie, xib, logr + +class compEbTreecorr: + + def __init__(self, x, y, dx, dy, rmin=5.0 / 3600.0, rmax=1.5, dlogr=0.05): + + self.catalog = treecorr.Catalog(x=x, y=y, v1=dx, v2=dy) + self.vvCorr = treecorr.VVCorrelation(min_sep=rmin, max_sep=rmax, bin_size=dlogr) + + def vcorr(self): + + self.vvCorr.process(self.catalog) + + def xiB(self, logr, xiplus, ximinus): + """ + Return estimate of pure B-mode correlation function + """ + # Integral of d(log r) ximinus(r) from r to infty: + dlogr = np.zeros_like(logr) + dlogr[1:-1] = 0.5 * (logr[2:] - logr[:-2]) + tmp = np.array(ximinus) * dlogr + integral = np.cumsum(tmp[::-1])[::-1] + return 0.5 * (xiplus - ximinus) + integral + + + def comp_eb(self): + self.vcorr() + xib = self.xiB(self.vvCorr.logr, self.vvCorr.xip, self.vvCorr.xim) + xie = self.vvCorr.xip - xib + return xie, xib, self.vvCorr.logr + + +def comp_eb_treecorr(u, v, du, dv, **kwargs): + """ + Compute E/B decomposition of astrometric error correlation function + + Parameters + ---------- + u, v : array_like. positions of objects. + du, dv : array_like. astrometric shift. + + returns + ------- + xie, xib, logr : array_like. E-mode, B-mode, + and log of binned distance separation in 2-point correlation function. + """ + cebt = compEbTreecorr(u, v, du, dv, **kwargs) + xie, xib, logr = cebt.comp_eb() + return xie, xib, logr From 51ed81b4794c772c277909c226cc347644685f97 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Leget Date: Wed, 26 Feb 2025 12:48:08 -0500 Subject: [PATCH 2/5] run black --- treegp/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/treegp/utils.py b/treegp/utils.py index ead06db..110caee 100644 --- a/treegp/utils.py +++ b/treegp/utils.py @@ -108,6 +108,7 @@ def comp_eb(u, v, du, dv, **kwargs): xie = xiplus - xib return xie, xib, logr + class compEbTreecorr: def __init__(self, x, y, dx, dy, rmin=5.0 / 3600.0, rmax=1.5, dlogr=0.05): @@ -130,7 +131,6 @@ def xiB(self, logr, xiplus, ximinus): integral = np.cumsum(tmp[::-1])[::-1] return 0.5 * (xiplus - ximinus) + integral - def comp_eb(self): self.vcorr() xib = self.xiB(self.vvCorr.logr, self.vvCorr.xip, self.vvCorr.xim) From 22395d12465d4971f1811e2824c9ced004d23666 Mon Sep 17 00:00:00 2001 From: PFLeget Date: Wed, 26 Feb 2025 12:51:17 -0500 Subject: [PATCH 3/5] Update _version.py --- treegp/_version.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/treegp/_version.py b/treegp/_version.py index 1e5422b..39424a3 100644 --- a/treegp/_version.py +++ b/treegp/_version.py @@ -1,2 +1,2 @@ -__version__ = "1.3.0" +__version__ = "1.2.0" __version_info__ = tuple(map(int, __version__.split("."))) From 23bcfd24e635f84e151ecbecbdd678b4cda3b8ff Mon Sep 17 00:00:00 2001 From: PFLeget Date: Wed, 26 Feb 2025 12:52:00 -0500 Subject: [PATCH 4/5] Update utils.py --- treegp/utils.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/treegp/utils.py b/treegp/utils.py index 110caee..1253832 100644 --- a/treegp/utils.py +++ b/treegp/utils.py @@ -29,13 +29,11 @@ def vcorr(x, y, dx, dy, rmin=5.0 / 3600.0, rmax=1.5, dlogr=0.05, maxpts=30000): if len(x) > maxpts: # Subsample array to get desired number of points rate = float(maxpts) / len(x) - print("Subsampling rate {:5.3f}%".format(rate * 100.0)) use = np.random.random(len(x)) <= rate x = x[use] y = y[use] dx = dx[use] dy = dy[use] - # print("Length ", len(x)) # Get index arrays that make all unique pairs i1, i2 = np.triu_indices(len(x)) # Omit self-pairs From 1897a7d0edc178ebce4a44b1e0552cd0450b6c60 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Leget Date: Wed, 26 Feb 2025 12:54:38 -0500 Subject: [PATCH 5/5] add treecorr eb comp test --- tests/test_hyp_search.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/test_hyp_search.py b/tests/test_hyp_search.py index 0419258..e3e7e86 100644 --- a/tests/test_hyp_search.py +++ b/tests/test_hyp_search.py @@ -134,6 +134,7 @@ def test_hyperparameter_search_2d(): if opt == "anisotropic": try: e_mode, b_mode, logr = treegp.comp_eb(x[:, 0], x[:, 1], y, y) + e_mode, b_mode, logr = treegp.comp_eb_treecorr(x[:, 0], x[:, 1], y, y) except: raise ValueError("Failed to compute E/B decomposition") # test if found hyperparameters are close the true hyperparameters.