From 78dcd4a9a2b5b0d872c231906f05e14765dede14 Mon Sep 17 00:00:00 2001 From: "Mathew S. Madhavacheril" Date: Mon, 5 Feb 2024 19:26:56 -0500 Subject: [PATCH] Update to v1.2 (#16) * delete data directory * allow arbitrary versions; change default to v1.2 * ACT DR4 calibration options * update to v1.2 * README note * make test use lensed Cls instead of unlensed ; update test values to reflect norm correction fix * increase tolerance for one test that seems to be slightly sensitive to CAMB version * version update --- .github/workflows/testing.yml | 4 +- README.md | 7 +- act_dr6_lenslike/__init__.py | 2 +- act_dr6_lenslike/act_dr6_lenslike.py | 104 ++++++++++++++++++++------ act_dr6_lenslike/data/README | 1 - act_dr6_lenslike/tests/test_act.py | 27 ++++--- act_dr6_lenslike/tests/test_cobaya.py | 38 +++------- get-act-data.sh | 10 +-- 8 files changed, 123 insertions(+), 70 deletions(-) delete mode 100644 act_dr6_lenslike/data/README diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml index 609dd90..d256870 100644 --- a/.github/workflows/testing.yml +++ b/.github/workflows/testing.yml @@ -35,8 +35,8 @@ jobs: uses: actions/cache@v3 id: cache-act with: - path: act_dr6_lenslike/data/v1.1 - key: ${{ runner.os }}-act-dr6-v1.1 + path: act_dr6_lenslike/data/v1.2 + key: ${{ runner.os }}-act-dr6-v1.2 - name: Download ACT DR6 Lensing Data if: steps.cache-act.outputs.cache-hit != 'true' diff --git a/README.md b/README.md index ce67c27..386cfe3 100644 --- a/README.md +++ b/README.md @@ -34,7 +34,7 @@ This can be performed automatically with the supplied `get-act-data.sh` script. Download the likelihood data tarball for ACT DR6 lensing from [NASA's LAMBDA archive](https://lambda.gsfc.nasa.gov/product/act/actadv_prod_table.html). -Extract the tarball into the `act_dr6_lenslike/data/` directory in the cloned repository such the directory `v1.1` is directly inside it. Only then should you proceed with the next steps. +Extract the tarball into the `act_dr6_lenslike/data/` directory in the cloned repository such the directory `v1.2` is directly inside it. Only then should you proceed with the next steps. ## Step 3: use in Python codes @@ -75,7 +75,10 @@ likelihood: ``` No other parameters need to be set. (e.g. do not manually set `like_corrections` or `no_like_corrections` here). -An example is provided in `ACTDR6LensLike-example.yaml` +An example is provided in `ACTDR6LensLike-example.yaml`. If, however, you are combining with +the ACT DR4 CMB 2-point power spectrum likelihood, you should also set `no_actlike_cmb_corrections: True` +(in addition to `lens_only: True` as described below). You do not need to do this if you are combining +with Planck CMB 2-point power spectrum likelihoods. ### Important parameters diff --git a/act_dr6_lenslike/__init__.py b/act_dr6_lenslike/__init__.py index a1c3ea9..e789285 100644 --- a/act_dr6_lenslike/__init__.py +++ b/act_dr6_lenslike/__init__.py @@ -1,5 +1,5 @@ """Likelihood for the Atacama Cosmology Telescope DR6 CMB lensing data.""" __author__ = ["Mathew Madhavacheril", "Ian Harrison"] -__version__ = "1.0.6" +__version__ = "1.2.0" from .act_dr6_lenslike import * diff --git a/act_dr6_lenslike/act_dr6_lenslike.py b/act_dr6_lenslike/act_dr6_lenslike.py index dcfac53..fe2775d 100644 --- a/act_dr6_lenslike/act_dr6_lenslike.py +++ b/act_dr6_lenslike/act_dr6_lenslike.py @@ -6,8 +6,7 @@ except: InstallableLikelihood = object import os -file_dir = os.path.abspath(os.path.dirname(__file__)) -data_dir = f"{file_dir}/data/v1.1/" +default_version = "v1.2" variants =[x.strip() for x in ''' act_baseline, @@ -51,7 +50,13 @@ def download(url, filename): return path def get_data(data_url="https://lambda.gsfc.nasa.gov/data/suborbital/ACT/ACT_dr6/likelihood/data/", - data_filename="ACT_dr6_likelihood_v1.1.tgz"): + data_filename_root="ACT_dr6_likelihood",version=None): + + if version is None: + version = default_version + data_filename = f"f{data_filename_root}_{version}.tgz" + file_dir = os.path.abspath(os.path.dirname(__file__)) + data_dir = f"{file_dir}/data/{version}/" if os.path.exists(os.path.join(file_dir, data_dir)): print('Data already exists at {}, not downloading again.'.format(os.path.join(file_dir, data_dir))) @@ -66,7 +71,7 @@ def get_data(data_url="https://lambda.gsfc.nasa.gov/data/suborbital/ACT/ACT_dr6/ download(data_url+data_filename, data_filename) tar = tarfile.open(data_filename) - tar.extractall(path=os.path.join(file_dir, data_dir).rstrip('v1.1/')) # this is not great + tar.extractall(path=os.path.join(file_dir, data_dir).rstrip(f'{version}/')) # this is not great tar.close() os.remove(data_filename) @@ -75,22 +80,45 @@ def get_data(data_url="https://lambda.gsfc.nasa.gov/data/suborbital/ACT/ACT_dr6/ def pp_to_kk(clpp,ell): return clpp * (ell*(ell+1.))**2. / 4. -def get_corrected_clkk(data_dict,clkk,cltt,clte,clee,clbb,suff=''): +def get_corrected_clkk(data_dict,clkk,cltt,clte,clee,clbb,suff='', + do_norm_corr=True, do_N1kk_corr=True, do_N1cmb_corr=True, + act_calib=False, no_like_cmb_corrections=False): + if no_like_cmb_corrections: + do_norm_corr = False + do_N1cmb_corr = False clkk_fid = data_dict['fiducial_cl_kk'] cl_dict = {'tt':cltt,'te':clte,'ee':clee,'bb':clbb} - N1_kk_corr = data_dict[f'dN1_kk{suff}'] @ (clkk-clkk_fid) + if do_N1kk_corr: + N1_kk_corr = data_dict[f'dN1_kk{suff}'] @ (clkk-clkk_fid) + else: + N1_kk_corr = 0 dNorm = data_dict[f'dAL_dC{suff}'] fid_norm = data_dict[f'fAL{suff}'] N1_cmb_corr = 0. norm_corr = 0. + + if act_calib and not('planck' in suff): + ocl = cl_dict['tt'] + fcl = data_dict[f'fiducial_cl_tt'] + ols = np.arange(ocl.size) + cal_ell_min = 1000 + cal_ell_max = 2000 + sel = np.s_[np.logical_and(ols>cal_ell_min,ols=2] = c[ls>=2] / fid_norm[ls>=2] - norm_corr = norm_corr + c + icl = cl_dict[s] + cldiff = ((icl/cal_fact)-data_dict[f'fiducial_cl_{s}']) + if do_N1cmb_corr: + N1_cmb_corr = N1_cmb_corr + (data_dict[f'dN1_{s}{suff}']@cldiff) + if do_norm_corr: + c = - 2. * (dNorm[i] @ cldiff) + if i==0: + ls = np.arange(c.size) + c[ls>=2] = c[ls>=2] / fid_norm[ls>=2] + norm_corr = norm_corr + c nclkk = clkk + norm_corr*clkk_fid + N1_kk_corr + N1_cmb_corr return nclkk @@ -200,10 +228,11 @@ def parse_variant(variant): chi_square = -2 lnlike """ -def load_data(variant, ddir=data_dir, +def load_data(variant, ddir=None, lens_only=False, apply_hartlap=True,like_corrections=True,mock=False, - nsims_act=796,nsims_planck=400,trim_lmax=2998,scale_cov=None): + nsims_act=796,nsims_planck=400,trim_lmax=2998,scale_cov=None, + version=None, act_cmb_rescale=False, act_calib=False): """ Given a data directory path, this function loads into a dictionary the data products necessary for evaluating the DR6 lensing likelihood. @@ -222,6 +251,12 @@ def load_data(variant, ddir=data_dir, """ # TODO: review defaults + if version is None: + version = default_version + + if ddir is None: + file_dir = os.path.abspath(os.path.dirname(__file__)) + ddir = f"{file_dir}/data/{version}/" if not os.path.exists(ddir): raise FileNotFoundError("Requested data directory {} does not exist.\ @@ -230,7 +265,9 @@ def load_data(variant, ddir=data_dir, with the act_dr6_lenslike.get_data() function.".format(ddir)) + print(f"Loading ACT DR6 lensing likelihood {version}...") v,baseline,include_planck = parse_variant(variant) + if include_planck and act_cmb_rescale: raise ValueError # output data @@ -281,6 +318,7 @@ def load_data(variant, ddir=data_dir, data_act = y[start:end].copy() d['data_binned_clkk'] = data_act nbins_act = data_act.size + binmat = np.loadtxt(f'{ddir}/binning_matrix_act.txt') d['full_binmat_act'] = binmat.copy() @@ -290,7 +328,22 @@ def load_data(variant, ddir=data_dir, d['binmat_act'] = standardize(ls,binmat[start:end,:],trim_lmax,extra_dims="xy") d['bcents_act'] = bcents[start:end].copy() + if act_cmb_rescale: + # load A_L_fid / A_L_ACT and standardize it + r = np.loadtxt("ratio_fid_over_act_wmap.txt") + rls = np.arange(r.size) + r[rls<2] = 0 + rs = standardize(rls,r,trim_lmax) + # bin it + rb = d['binmat_act'] @ rs + # correct data + print("Binned ACT rescaling corrections: ", rb) + d['data_binned_clkk'] = d['data_binned_clkk'] / rb**2 + + if lens_only: + if act_cmb_rescale: raise ValueError + if act_calib: raise ValueError if include_planck: if v not in [None,'cinpaint']: raise ValueError(f"Combination of {v} with Planck is not available") fcov = np.loadtxt(f'{ddir}/covmat_actplanck_cmbmarg.txt') @@ -360,8 +413,7 @@ def load_data(variant, ddir=data_dir, nsims = min(nsims_act,nsims_planck) if include_planck else nsims_act hartlap_correction = (nsims-nbins-2.)/(nsims-1.) if apply_hartlap: - pass - #warnings.warn("Hartlap correction to cinv: ", hartlap_correction) + warnings.warn(f"Hartlap correction to cinv: {hartlap_correction}") else: warnings.warn(f"Disabled Hartlap correction to cinv: {hartlap_correction}") hartlap_correction = 1.0 @@ -381,7 +433,8 @@ def load_data(variant, ddir=data_dir, return d -def generic_lnlike(data_dict,ell_kk,cl_kk,ell_cmb,cl_tt,cl_ee,cl_te,cl_bb,trim_lmax=2998,return_theory=False): +def generic_lnlike(data_dict,ell_kk,cl_kk,ell_cmb,cl_tt,cl_ee,cl_te,cl_bb,trim_lmax=2998, + return_theory=False,do_norm_corr=True,act_calib=False,no_actlike_cmb_corrections=False): cl_kk = standardize(ell_kk,cl_kk,trim_lmax) cl_tt = standardize(ell_cmb,cl_tt,trim_lmax) @@ -391,7 +444,9 @@ def generic_lnlike(data_dict,ell_kk,cl_kk,ell_cmb,cl_tt,cl_ee,cl_te,cl_bb,trim_l d = data_dict cinv = d['cinv'] - clkk_act = get_corrected_clkk(data_dict,cl_kk,cl_tt,cl_te,cl_ee,cl_bb) if d['likelihood_corrections'] else cl_kk + clkk_act = get_corrected_clkk(data_dict,cl_kk,cl_tt,cl_te,cl_ee,cl_bb, + do_norm_corr=do_norm_corr,act_calib=act_calib, + no_like_cmb_corrections=no_actlike_cmb_corrections) if d['likelihood_corrections'] else cl_kk bclkk = d['binmat_act'] @ clkk_act if d['include_planck']: clkk_planck = get_corrected_clkk(data_dict,cl_kk,cl_tt,cl_te,cl_ee,cl_bb,'_planck') if d['likelihood_corrections'] else cl_kk @@ -416,6 +471,7 @@ class ACTDR6LensLike(InstallableLikelihood): nsims_act = 792. # Number of sims used for covmat; used in Hartlap correction nsims_planck = 400. # Number of sims used for covmat; used in Hartlap correction no_like_corrections = False + no_actlike_cmb_corrections = False lens_only = False # Any ells above this will be discarded; likelihood must at least request ells up to this trim_lmax = 2998 @@ -428,6 +484,9 @@ class ACTDR6LensLike(InstallableLikelihood): zmax = None scale_cov = None varying_cmb_alens = False # Whether to divide the theory spectrum by Alens + version = None + act_cmb_rescale = False + act_calib = False def initialize(self): if self.lens_only: self.no_like_corrections = True @@ -435,7 +494,8 @@ def initialize(self): self.data = load_data(variant=self.variant,lens_only=self.lens_only, like_corrections=not(self.no_like_corrections),apply_hartlap=self.apply_hartlap, mock=self.mock,nsims_act=self.nsims_act,nsims_planck=self.nsims_planck, - trim_lmax=self.trim_lmax,scale_cov=self.scale_cov) + trim_lmax=self.trim_lmax,scale_cov=self.scale_cov,version=self.version, + act_cmb_rescale=self.act_cmb_rescale,act_calib=self.act_calib) if self.no_like_corrections: self.requested_cls = ["pp"] @@ -474,7 +534,9 @@ def loglike(self, cl, **params_values): else: cl_kk = pp_to_kk(clpp,ell) - logp = generic_lnlike(self.data,ell,cl_kk,ell,cl['tt'],cl['ee'],cl['te'],cl['bb'],self.trim_lmax) + logp = generic_lnlike(self.data,ell,cl_kk,ell,cl['tt'],cl['ee'],cl['te'],cl['bb'],self.trim_lmax, + do_norm_corr=not(self.act_cmb_rescale),act_calib=self.act_calib, + no_actlike_cmb_corrections=self.no_actlike_cmb_corrections) self.log.debug( f"ACT-DR6-lensing-like lnLike value = {logp} (chisquare = {-2 * logp})") return logp diff --git a/act_dr6_lenslike/data/README b/act_dr6_lenslike/data/README deleted file mode 100644 index e72461d..0000000 --- a/act_dr6_lenslike/data/README +++ /dev/null @@ -1 +0,0 @@ -Directory for data files \ No newline at end of file diff --git a/act_dr6_lenslike/tests/test_act.py b/act_dr6_lenslike/tests/test_act.py index 96cd86c..8e6a2b3 100644 --- a/act_dr6_lenslike/tests/test_act.py +++ b/act_dr6_lenslike/tests/test_act.py @@ -1,30 +1,33 @@ import unittest import act_dr6_lenslike as alike import numpy as np -data_dir = alike.data_dir +import os +file_dir = os.path.abspath(os.path.dirname(__file__)) +version = alike.default_version +data_dir = f"{file_dir}/../data/{version}/" class ACTLikeTest(unittest.TestCase): def generic_call(self,variant,lens_only,exp_chisq=None,return_theory=False): - data_file = data_dir+'like_corrs/cosmo2017_10K_acc3_lenspotentialCls.dat' try: - ell, cl_tt, cl_ee, cl_bb, cl_te, cl_pp, cl_tp, cl_ep= np.loadtxt(data_file, unpack=True) + ell, cl_tt, cl_ee, cl_bb, cl_te = np.loadtxt(data_dir+'like_corrs/cosmo2017_10K_acc3_lensedCls.dat', unpack=True) + ellp, _, _, _, _, cl_pp, _, _= np.loadtxt(data_dir+'like_corrs/cosmo2017_10K_acc3_lenspotentialCls.dat', unpack=True) except OSError: - raise - finally: print('Required data file not found at {}'.format(data_file)) print('Please obtain it and place it correctly.') print('The script get-act-data.sh will download and place it.') + raise + prefac = 2*np.pi/ell/(ell+1.) cl_kk=cl_pp/4*2*np.pi cl_bb = cl_bb*prefac cl_tt = cl_tt*prefac cl_ee = cl_ee*prefac cl_te = cl_te*prefac - data_dict = alike.load_data(variant,lens_only=lens_only,like_corrections=not(lens_only)) - ell_kk = ell + data_dict = alike.load_data(variant,lens_only=lens_only,like_corrections=not(lens_only),version=version) + ell_kk = ellp ell_cmb=ell if return_theory: @@ -37,21 +40,21 @@ def generic_call(self,variant,lens_only,exp_chisq=None,return_theory=False): def test_act_baseline_lensonly(self): self.generic_call('act_baseline',True,14.06) def test_act_baseline(self): - self.generic_call('act_baseline',False,14.71) + self.generic_call('act_baseline',False,14.13) def test_act_baseline_return_theory(self): - self.generic_call('act_baseline',False,14.71,return_theory=True) + self.generic_call('act_baseline',False,14.13,return_theory=True) def test_actplanck_baseline_lensonly(self): self.generic_call('actplanck_baseline',True,21.07) def test_actplanck_baseline(self): - self.generic_call('actplanck_baseline',False,21.97) + self.generic_call('actplanck_baseline',False,21.46) def test_act_extended_lensonly(self): self.generic_call('act_extended',True,17.66) def test_act_extended(self): - self.generic_call('act_extended',False,17.26) + self.generic_call('act_extended',False,17.72) def test_actplanck_extended_lensonly(self): self.generic_call('actplanck_extended',True,24.4) def test_actplanck_extended(self): - self.generic_call('actplanck_extended',False,24.26) + self.generic_call('actplanck_extended',False,24.75) def test_act_polonly_lensonly(self): self.generic_call('act_polonly',True,309.71) def test_act_cibdeproj_lensonly(self): diff --git a/act_dr6_lenslike/tests/test_cobaya.py b/act_dr6_lenslike/tests/test_cobaya.py index e0072c0..6d07359 100644 --- a/act_dr6_lenslike/tests/test_cobaya.py +++ b/act_dr6_lenslike/tests/test_cobaya.py @@ -16,7 +16,9 @@ 'omnuh2': 0.00064, }, "theory" : { - "camb" : { 'extra_args' : { + "camb" : { + 'ignore_obsolete': True, + 'extra_args' : { 'lmax' : 10000, 'lens_margin' : 1250, 'lens_potential_accuracy' : 4, @@ -31,21 +33,11 @@ class ACTLikeTest(unittest.TestCase): - # Note that tests marked: 'fail on new camb' have had the chi2 revised - # from the original number which was with respect to a simulated spectrum - # from this analysis: https://mapsims.readthedocs.io/en/latest/camb.html - # This spectrum is included as like_corrs/cosmo2017_10K_acc3_lenspotentialCls.dat - # and is used for the tests for the base python likelihood. - # If anyone can reproduce these numbers with settings for a modern camb, - # please get in touch! - # We do not believe this affects the validity of the results, as the chi2 - # differences are small, and the previously compared spectrum was somewhat - # arbitray. def initialize(self): install({"likelihood": {"act_dr6_lenslike.ACTDR6LensLike": None}}) - def generic_call(self,variant,lens_only,exp_chisq=None): + def generic_call(self,variant,lens_only,exp_chisq=None,places=1): info['likelihood'] = {'ACTDR6LensLike' : {'external' : ACTDR6LensLike, 'variant' : variant, @@ -56,44 +48,38 @@ def generic_call(self,variant,lens_only,exp_chisq=None): chisq = -2 * loglikes[0] - self.assertAlmostEqual(chisq, exp_chisq, 1) + self.assertAlmostEqual(chisq, exp_chisq, places) def test_act_baseline_lensonly(self): - self.generic_call('act_baseline',True,14.06) + self.generic_call('act_baseline',True,14.00) def test_act_baseline(self): - # self.generic_call('act_baseline',False,14.71) # fail on new camb - self.generic_call('act_baseline',False,14.07) + self.generic_call('act_baseline',False,14.00) def test_actplanck_baseline_lensonly(self): - # self.generic_call('actplanck_baseline',True,21.07) # fail on new camb self.generic_call('actplanck_baseline',True,20.97) def test_actplanck_baseline(self): - # self.generic_call('actplanck_baseline',False,21.97) # fail on new camb - self.generic_call('actplanck_baseline',False,21.36) + self.generic_call('actplanck_baseline',False,21.20) def test_act_extended_lensonly(self): - # self.generic_call('act_extended',True,17.66) # fail on new camb self.generic_call('act_extended',True,17.84) def test_act_extended(self): - # self.generic_call('act_extended',False,17.26) # fail on new camb - self.generic_call('act_extended',False,17.94) + # This seems to have extra sensitivity to CAMB version + self.generic_call('act_extended',False,17.94,places=0) def test_actplanck_extended_lensonly(self): - # self.generic_call('actplanck_extended',True,24.4) # fail on new camb self.generic_call('actplanck_extended',True,24.52) def test_actplanck_extended(self): - # self.generic_call('actplanck_extended',False,24.26) # fail on new camb - self.generic_call('actplanck_extended',False,24.91) + self.generic_call('actplanck_extended',False,24.77) def test_act_polonly_lensonly(self): self.generic_call('act_polonly',True,309.71) def test_act_cibdeproj_lensonly(self): - self.generic_call('act_cibdeproj',True,15.16) + self.generic_call('act_cibdeproj',True,15.10) def test_act_cinpaint_lensonly(self): self.generic_call('act_cinpaint',True,15.94) diff --git a/get-act-data.sh b/get-act-data.sh index d398386..08adc36 100755 --- a/get-act-data.sh +++ b/get-act-data.sh @@ -2,7 +2,7 @@ # script from Joe Zuntz -if [ -d act_dr6_lenslike/data/v1.1 ] +if [ -d act_dr6_lenslike/data/v1.2 ] then echo ACT DR6 Lensing data already downloaded elif ! command -v wget &> /dev/null @@ -11,8 +11,8 @@ then else mkdir -p act_dr6_lenslike/data pushd act_dr6_lenslike/data - wget https://lambda.gsfc.nasa.gov/data/suborbital/ACT/ACT_dr6/likelihood/data/ACT_dr6_likelihood_v1.1.tgz - tar -zxvf ACT_dr6_likelihood_v1.1.tgz - rm ACT_dr6_likelihood_v1.1.tgz + wget https://lambda.gsfc.nasa.gov/data/suborbital/ACT/ACT_dr6/likelihood/data/ACT_dr6_likelihood_v1.2.tgz + tar -zxvf ACT_dr6_likelihood_v1.2.tgz + rm ACT_dr6_likelihood_v1.2.tgz popd -fi \ No newline at end of file +fi