From 94e7b75acc7262891a953ee526bcd6a347beb8fb Mon Sep 17 00:00:00 2001 From: Joe Zuntz Date: Thu, 11 May 2023 16:39:03 +0100 Subject: [PATCH 01/12] Add ACT DR6 likelihood, first attempt --- examples/act-dr6-lens-priors.ini | 3 + examples/act-dr6-lens-values.ini | 24 + examples/act-dr6-lens.ini | 45 ++ likelihood/act-dr6-lens/.gitignore | 1 + likelihood/act-dr6-lens/LICENSE | 24 + likelihood/act-dr6-lens/act_dr6_lenslike.py | 490 ++++++++++++++++++ .../act_dr6_lenslike_interface.py | 42 ++ likelihood/act-dr6-lens/get-act-data.sh | 14 + 8 files changed, 643 insertions(+) create mode 100644 examples/act-dr6-lens-priors.ini create mode 100644 examples/act-dr6-lens-values.ini create mode 100644 examples/act-dr6-lens.ini create mode 100644 likelihood/act-dr6-lens/.gitignore create mode 100644 likelihood/act-dr6-lens/LICENSE create mode 100644 likelihood/act-dr6-lens/act_dr6_lenslike.py create mode 100644 likelihood/act-dr6-lens/act_dr6_lenslike_interface.py create mode 100755 likelihood/act-dr6-lens/get-act-data.sh diff --git a/examples/act-dr6-lens-priors.ini b/examples/act-dr6-lens-priors.ini new file mode 100644 index 00000000..f0921feb --- /dev/null +++ b/examples/act-dr6-lens-priors.ini @@ -0,0 +1,3 @@ +[cosmological_parameters] +ombh2 = normal 0.02233 0.00036 +n_s = normal 0.96 0.02 \ No newline at end of file diff --git a/examples/act-dr6-lens-values.ini b/examples/act-dr6-lens-values.ini new file mode 100644 index 00000000..4af5382a --- /dev/null +++ b/examples/act-dr6-lens-values.ini @@ -0,0 +1,24 @@ + +[cosmological_parameters] +omch2 = 0.005 0.11 0.99 +ombh2 = 0.01 0.03 0.02233 +log1e10As = 1.61 3.0448 4.0 +n_s = 0.86 0.96 1.06 +cosmomc_theta = 0.005 0.0104 0.1 + +omega_k = 0.0 ;spatial curvature + +;neutrinos +mnu = 0.06 +nnu = 3.046 +num_massive_neutrinos = 1 + +;helium +yhe = 0.245341 ;helium mass fraction + +;reionization +tau = 0.0543 ;reionization optical depth + +[halo_model_parameters] +A = 3.13 +eta = 0.603 \ No newline at end of file diff --git a/examples/act-dr6-lens.ini b/examples/act-dr6-lens.ini new file mode 100644 index 00000000..cc7135b1 --- /dev/null +++ b/examples/act-dr6-lens.ini @@ -0,0 +1,45 @@ +[runtime] +sampler = test +root = ${PWD} + +[test] +save_dir=output/act_dr6_lens +fatal_errors=T + + +[pipeline] +; these names refer to sections later in the file: +modules = consistency camb act_dr6_lens +values = examples/act-dr6-lens-values.ini +priors = examples/act-dr6-lens-priors.ini +quiet=F +debug=T +timing=F + + +[act_dr6_lens] +file = ./likelihood/act-dr6-lens/act_dr6_lenslike_interface.py + + + +; The consistency module translates between our chosen parameterization +; and any other that modules in the pipeline may want (e.g. camb) +[consistency] +file = ./utility/consistency/consistency_interface.py +cosmomc_theta = T + +[camb] +file = boltzmann/camb/camb_interface.py +mode = cmb +lmax = 4000 ;max ell to use for cmb calculation +lens_margin = 1250 +lens_potential_accuracy = 4 +feedback=2 ;amount of output to print +AccuracyBoost=1.0 ;CAMB accuracy boost parameter +lSampleBoost = 1.0 +lAccuracyBoost = 1.0 +do_tensors = True ;include tensor modes +do_lensing = true ;lensing is required w/ Planck data +NonLinear = lens +theta_H0_range = "20 100" +halofit_version = mead2016 diff --git a/likelihood/act-dr6-lens/.gitignore b/likelihood/act-dr6-lens/.gitignore new file mode 100644 index 00000000..8fce6030 --- /dev/null +++ b/likelihood/act-dr6-lens/.gitignore @@ -0,0 +1 @@ +data/ diff --git a/likelihood/act-dr6-lens/LICENSE b/likelihood/act-dr6-lens/LICENSE new file mode 100644 index 00000000..5508fc2c --- /dev/null +++ b/likelihood/act-dr6-lens/LICENSE @@ -0,0 +1,24 @@ +BSD 2-Clause License + +Copyright (c) 2023, Members of the Atacama Cosmology Telescope Collaboration + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + +1. Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + +2. Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE +FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/likelihood/act-dr6-lens/act_dr6_lenslike.py b/likelihood/act-dr6-lens/act_dr6_lenslike.py new file mode 100644 index 00000000..222c69e0 --- /dev/null +++ b/likelihood/act-dr6-lens/act_dr6_lenslike.py @@ -0,0 +1,490 @@ +__author__ = "Mathew Madhavacheril" +__version__ = "1.0.0" + +import numpy as np +import warnings +from scipy.interpolate import interp1d +InstallableLikelihood = object +import os +file_dir = os.path.abspath(os.path.dirname(__file__)) +data_dir = f"{file_dir}/data/v1.1/" + +variants =[x.strip() for x in ''' +act_baseline, +act_extended, +actplanck_baseline, +actplanck_extended, +act_polonly, +act_cibdeproj, +act_cinpaint +'''.strip().replace('\n','').split(',')] + + +# ================ +# HELPER FUNCTIONS +# ================ + +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=''): + 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) + dNorm = data_dict[f'dAL_dC{suff}'] + fid_norm = data_dict[f'fAL{suff}'] + N1_cmb_corr = 0. + norm_corr = 0. + for i,s in enumerate(['tt','ee','bb','te']): + cldiff = (cl_dict[s]-data_dict[f'fiducial_cl_{s}']) + N1_cmb_corr = N1_cmb_corr + (data_dict[f'dN1_{s}{suff}']@cldiff) + 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 + +def standardize(ls,cls,trim_lmax,lbuffer=2,extra_dims="y"): + cstart = int(ls[0]) + diffs = np.diff(ls) + if not(np.all(np.isclose(diffs,1.))): raise ValueError("Multipoles are not spaced by 1") + if not(cstart<=2): raise ValueError("Multipoles start at value greater than 2") + nlen = trim_lmax+lbuffer + cend = nlen - cstart + if extra_dims=="xyy": + out = np.zeros((cls.shape[0],nlen,nlen)) + out[:,cstart:,cstart:] = cls[:,:cend,:cend] + elif extra_dims=="yy": + out = np.zeros((nlen,nlen)) + out[cstart:,cstart:] = cls[:cend,:cend] + elif extra_dims=="xy": + out = np.zeros((cls.shape[0],nlen)) + out[:,cstart:] = cls[:,:cend] + elif extra_dims=="y": + out = np.zeros(nlen) + out[cstart:] = cls[:cend] + else: + raise ValueError + return out + +def get_limber_clkk_flat_universe(results,Pfunc,lmax,kmax,nz,zsrc=None): + # Adapting code from Antony Lewis' CAMB notebook + if zsrc is None: + chistar = results.conformal_time(0)- results.tau_maxvis + else: + chistar = results.comoving_radial_distance(zsrc) + chis = np.linspace(0,chistar,nz) + zs=results.redshift_at_comoving_radial_distance(chis) + dchis = (chis[2:]-chis[:-2])/2 + chis = chis[1:-1] + zs = zs[1:-1] + + #Get lensing window function (flat universe) + win = ((chistar-chis)/(chis**2*chistar))**2 + #Do integral over chi + ls = np.arange(0,lmax+2, dtype=np.float64) + cl_kappa=np.zeros(ls.shape) + w = np.ones(chis.shape) #this is just used to set to zero k values out of range of interpolation + for i, l in enumerate(ls[2:]): + k=(l+0.5)/chis + w[:]=1 + w[k<1e-4]=0 + w[k>=kmax]=0 + cl_kappa[i+2] = np.dot(dchis, w*Pfunc.P(zs, k, grid=False)*win/k**4) + cl_kappa*= (ls*(ls+1))**2 + return cl_kappa + +def get_camb_lens_obj(nz,kmax,zmax=None): + import camb + pars = camb.CAMBparams() + # This cosmology is purely to go from chis->zs for limber integration; + # the details do not matter + pars.set_cosmology(H0=67.5, ombh2=0.022, omch2=0.122) + pars.InitPower.set_params(ns=0.965) + results= camb.get_background(pars) + nz = nz + if zmax is None: + chistar = results.conformal_time(0)- results.tau_maxvis + else: + chistar = results.comoving_radial_distance(zmax) + chis = np.linspace(0,chistar,nz) + zs=results.redshift_at_comoving_radial_distance(chis) + cobj = {"CAMBdata": None, + "Pk_interpolator": { "z": zs, + "k_max": kmax, + "nonlinear": True, + "vars_pairs": ([["Weyl", "Weyl"]])}} + return cobj + + +def parse_variant(variant): + + variant = variant.lower().strip() + if variant not in variants: raise ValueError + + v = None + if '_extended' in variant: + baseline = False + else: + baseline = True + if '_baseline' not in variant: + v = variant.split('_')[-1] + + include_planck = True if 'actplanck' in variant else False + return v,baseline,include_planck + + +# ================== +# Generic likelihood +# ================== + +""" +data_dict = load_data(data_directory) # pre-load data +# for each predicted spectra in chain +# cl_kk is CMB lensing convergence power spectrum (dimensionless, +# no ell or 2pi factors) +# cl_tt, cl_ee, cl_te, cl_bb are lensed CMB power spectra +# (muK^2 units, no ell or 2pi factors) +lnlike = generic_lnlike(data_dict,cl_kk,cl_tt,cl_ee,cl_te,cl_bb) +This returns ln(Likelihood) +so for example, +chi_square = -2 lnlike +""" + +def load_data(variant,lens_only=False, + apply_hartlap=True,like_corrections=True,mock=False, + nsims_act=796,nsims_planck=400,trim_lmax=2998,scale_cov=None): + """ + Given a data directory path, this function loads into a dictionary + the data products necessary for evaluating the DR6 lensing likelihood. + This includes: + 1. the ACT lensing bandpowers. Planck lensing bandpowers will be + appended if include_planck is True. + 2. the associated binning matrix to be applied to a theory curve + 3. the associated covariance matrix + 4. data products associated with applying likelihood corrections + + All these products will be standardized so that they apply + to theory curves specified from L=0 to trim_lmax. + + A Hartlap correction will be applied to the covariance matrix + corresponding to the lower of the number of simulations involved. + + """ + # TODO: review defaults + + + v,baseline,include_planck = parse_variant(variant) + + + # output data + d = {} + + if lens_only and like_corrections: raise ValueError("Likelihood corrections should not be used in lens_only runs.") + if not(lens_only) and not(like_corrections): + warnings.warn("Neither using CMB-marginalized covariance matrix nor including likelihood corrections. Effective covariance may be underestimated.") + + d['include_planck'] = include_planck + d['likelihood_corrections'] = like_corrections + + # Fiducial spectra + if like_corrections: + f_ls, f_tt, f_ee, f_bb, f_te = np.loadtxt(f"{data_dir}/like_corrs/cosmo2017_10K_acc3_lensedCls.dat",unpack=True) + f_tt = f_tt / (f_ls * (f_ls+1.)) * 2. * np.pi + f_ee = f_ee / (f_ls * (f_ls+1.)) * 2. * np.pi + f_bb = f_bb / (f_ls * (f_ls+1.)) * 2. * np.pi + f_te = f_te / (f_ls * (f_ls+1.)) * 2. * np.pi + + fd_ls, f_dd = np.loadtxt(f"{data_dir}/like_corrs/cosmo2017_10K_acc3_lenspotentialCls.dat",unpack=True,usecols=[0,5]) + f_kk = f_dd * 2. * np.pi / 4. + d['fiducial_cl_tt'] = standardize(f_ls,f_tt,trim_lmax) + d['fiducial_cl_te'] = standardize(f_ls,f_te,trim_lmax) + d['fiducial_cl_ee'] = standardize(f_ls,f_ee,trim_lmax) + d['fiducial_cl_bb'] = standardize(f_ls,f_bb,trim_lmax) + d['fiducial_cl_kk'] = standardize(fd_ls,f_kk,trim_lmax) + + + # Return data bandpowers, covariance matrix and binning matrix + ddir = data_dir + if baseline: + start = 2 + end = -6 + else: + start = 2 + end = -3 + + if v is None: + y = np.loadtxt(f'{ddir}/clkk_bandpowers_act.txt') + elif v=='cinpaint': + y = np.loadtxt(f'{ddir}/clkk_bandpowers_act_cinpaint.txt') + elif v=='polonly': + y = np.loadtxt(f'{ddir}/clkk_bandpowers_act_polonly.txt') + elif v=='cibdeproj': + y = np.loadtxt(f'{ddir}/clkk_bandpowers_act_cibdeproj.txt') + nbins_tot_act = y.size + d['full_data_binned_clkk_act'] = y.copy() + 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() + pells = np.arange(binmat.shape[1]) + bcents = binmat@pells + ls = np.arange(binmat.shape[1]) + d['binmat_act'] = standardize(ls,binmat[start:end,:],trim_lmax,extra_dims="xy") + d['bcents_act'] = bcents.copy() + + if lens_only: + 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') + else: + if v=='cibdeproj': + fcov = np.loadtxt(f"{ddir}/covmat_act_cibdeproj_cmbmarg.txt") + elif v=='pol': + fcov = np.loadtxt(f"{ddir}/covmat_act_polonly_cmbmarg.txt") + else: + fcov = np.loadtxt(f"{ddir}/covmat_act_cmbmarg.txt") + else: + if v not in [None,'cinpaint']: raise ValueError(f"Covmat for {v} without CMB marginalization is not available") + if include_planck: + fcov = np.loadtxt(f'{ddir}/covmat_actplanck.txt') + else: + fcov = np.loadtxt(f'{ddir}/covmat_act.txt') + + d['full_act_cov'] = fcov.copy() + # Remove trailing bins from ACT part + sel = np.s_[nbins_tot_act+end:nbins_tot_act] + cov = np.delete(np.delete(fcov,sel,0),sel,1) + # Remove leading bins from ACT part + sel = np.s_[:start] + cov = np.delete(np.delete(cov,sel,0),sel,1) + + # Test + covmat = np.loadtxt(f'{ddir}/covmat_act.txt') + covmat1 = covmat[start:end,start:end] + cdiff = cov[:nbins_act,:nbins_act] - covmat1 + if not(np.all(np.isclose(cdiff,0))): raise ValueError + + if include_planck: + data_planck = np.loadtxt(f'{ddir}/clkk_bandpowers_planck.txt') + d['data_binned_clkk'] = np.append(d['data_binned_clkk'],data_planck) + binmat = np.loadtxt(f'{ddir}/binning_matrix_planck.txt') + pells = np.arange(binmat.shape[1]) + bcents = binmat@pells + ls = np.arange(binmat.shape[1]) + d['binmat_planck'] = standardize(ls,binmat,trim_lmax,extra_dims="xy") + d['bcents_planck'] = bcents.copy() + + if like_corrections: + # Load matrices + cmat = np.load(f"{ddir}/like_corrs/norm_correction_matrix_Lmin0_Lmax4000.npy") + ls = np.arange(cmat.shape[1]) + d['dAL_dC'] = standardize(ls,cmat,trim_lmax,extra_dims="xyy") + if include_planck: + cmat = np.load(f"{ddir}/like_corrs/P18_norm_correction_matrix_Lmin0_Lmax3000.npy") + ls = np.arange(cmat.shape[1]) + d['dAL_dC_planck'] = standardize(ls,cmat,trim_lmax,extra_dims="xyy") + + + fAL_ls,fAL = np.loadtxt(f"{ddir}/like_corrs/n0mv_fiducial_lmin600_lmax3000_Lmin0_Lmax4000.txt") + d['fAL'] = standardize(fAL_ls,fAL,trim_lmax,extra_dims="y") + if include_planck: + fAL_ls,fAL = np.loadtxt(f"{ddir}/like_corrs/PLANCK_n0mv_fiducial_lmin600_lmax3000_Lmin0_Lmax3000.txt") + d['fAL_planck'] = standardize(fAL_ls,fAL,trim_lmax,extra_dims="y") + + for spec in ['kk','tt','ee','bb','te']: + n1mat = np.loadtxt(f"{ddir}/like_corrs/N1der_{spec.upper()}_lmin600_lmax3000_full.txt") + d[f'dN1_{spec}'] = standardize(fAL_ls,n1mat,trim_lmax,extra_dims="yy") + if include_planck: + n1mat = np.loadtxt(f"{ddir}/like_corrs/N1_planck_der_{spec.upper()}_lmin100_lmax2048.txt") + d[f'dN1_{spec}_planck'] = standardize(fAL_ls,n1mat,trim_lmax,extra_dims="yy") + + nbins = d['data_binned_clkk'].size + 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) + else: + warnings.warn(f"Disabled Hartlap correction to cinv: {hartlap_correction}") + hartlap_correction = 1.0 + if scale_cov is not None: + warnings.warn(f"Covariance has been artificially scaled by: {scale_cov}") + cov = cov * scale_cov + d['cov'] = cov + cinv = np.linalg.inv(cov) * hartlap_correction + d['cinv'] = cinv + + if mock: + mclpp = np.loadtxt(f"{self.ddir}/cls_default_dr6_accuracy.txt",usecols=[5]) + ls = np.arange(mclpp.size) + mclkk = mclpp * 2. * np.pi / 4. + self.clkk_data = self.binning_matrix @ mclkk[:self.kLmax] + + return d + + +def generic_lnlike(data_dict,ell_kk,cl_kk,ell_cmb,cl_tt,cl_ee,cl_te,cl_bb,trim_lmax = 2998): + + cl_kk = standardize(ell_kk,cl_kk,trim_lmax) + cl_tt = standardize(ell_cmb,cl_tt,trim_lmax) + cl_ee = standardize(ell_cmb,cl_ee,trim_lmax) + cl_bb = standardize(ell_cmb,cl_bb,trim_lmax) + cl_te = standardize(ell_cmb,cl_te,trim_lmax) + + 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 + 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 + bclkk = np.append(bclkk, d['binmat_planck'] @ clkk_planck) + delta = d['data_binned_clkk'] - bclkk + lnlike = -0.5 * np.dot(delta,np.dot(cinv,delta)) + return lnlike + + +# ================= +# Cobaya likelihood +# ================= + +class GenericLimberCosmicShear(InstallableLikelihood): + + zsrc: float + ngal_arcmin2: float + fsky: float + glmin = 10 + lmin = 10 + lmax = 500 + nell = 20 + shape_std = 0.3 + draw_noise = False + nz = 40 + trim_lmax = 599 + kmax = 10 + cmb_noise = None + + def initialize(self): + from orphics import stats + import pyfisher + bin_edges = np.geomspace(self.glmin,self.lmax,self.nell) + bin_edges = bin_edges[bin_edges>self.lmin] + self.binner = stats.bin1D(bin_edges) + self.ls = np.arange(0,self.trim_lmax+2) + if self.cmb_noise is None: + nls_dict = {'kk': lambda x: x*0+self.shape_std**2/(2.*self.ngal_arcmin2*1.18e7)} + else: + import pyfisher + ls,nls = pyfisher.get_lensing_nl(self.cmb_noise) + nls_dict = {'kk': interp1d(ls,nls,bounds_error=True)} + cl_kk = self.get_mock_theory() + cls_dict = {'kk': interp1d(self.ls,cl_kk)} + self.d = {} + self.d['data_binned_clkk'] = self.binner.bin(self.ls,cl_kk)[1] + cov = pyfisher.gaussian_band_covariance(bin_edges,['kk'],cls_dict,nls_dict,interpolate=False)[:,0,0] / self.fsky + cinv = np.diag(1./cov) + self.d['cinv'] = cinv + + def get_mock_theory(self): + import camb + from camb import model + pars = camb.CAMBparams() + pars.set_cosmology(H0=67.5, ombh2=0.022, omch2=0.122) + pars.InitPower.set_params(ns=0.965) + results = camb.get_background(pars) + PK = camb.get_matter_power_interpolator(pars, nonlinear=True, + hubble_units=False, k_hunit=False, kmax=self.kmax, + var1=model.Transfer_Weyl,var2=model.Transfer_Weyl, zmax=self.zsrc) + return get_limber_clkk_flat_universe(results,PK,self.trim_lmax,self.kmax,self.nz,zsrc=self.zsrc) + + + def get_requirements(self): + cobj = get_camb_lens_obj(self.nz,self.kmax,self.zsrc) + return cobj + + def logp(self, **params_values): + cl_kk = self.get_limber_clkk( **params_values) + bclkk = self.binner.bin(self.ls,cl_kk)[1] + delta = self.d['data_binned_clkk'] - bclkk + logp = -0.5 * np.dot(delta,np.dot(self.d['cinv'],delta)) + self.log.debug( + f"Generic Cosmic Shear lnLike value = {logp} (chisquare = {-2 * logp})") + return logp + + def get_limber_clkk(self,**params_values): + Pfunc = self.theory.get_Pk_interpolator(var_pair=("Weyl", "Weyl"), nonlinear=True, extrap_kmax=30.) + results = self.provider.get_CAMBdata() + return get_limber_clkk_flat_universe(results,Pfunc,self.trim_lmax,self.kmax,self.nz,zsrc=self.zsrc) + +class ACTDR6LensLike(InstallableLikelihood): + + lmax: int + mock = False + 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 + lens_only = False + # Any ells above this will be discarded; likelihood must at least request ells up to this + trim_lmax = 2998 + variant = "act_baseline" + apply_hartlap = True + # Limber integral parameters + limber = False + nz = 100 + kmax = 10 + scale_cov = None + alens = False # Whether to divide the theory spectrum by Alens + + def initialize(self): + if self.lens_only: self.no_like_corrections = True + if self.lmax Date: Thu, 11 May 2023 16:45:30 +0100 Subject: [PATCH 02/12] fix popd arg --- likelihood/act-dr6-lens/get-act-data.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/likelihood/act-dr6-lens/get-act-data.sh b/likelihood/act-dr6-lens/get-act-data.sh index cc462264..243ce94c 100755 --- a/likelihood/act-dr6-lens/get-act-data.sh +++ b/likelihood/act-dr6-lens/get-act-data.sh @@ -10,5 +10,5 @@ else 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 - popd data + popd fi \ No newline at end of file From b9db0612a5210c579b137b984755bcbdd3364516 Mon Sep 17 00:00:00 2001 From: Joe Zuntz Date: Thu, 11 May 2023 16:47:58 +0100 Subject: [PATCH 03/12] fix prior --- examples/act-dr6-lens-values.ini | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/act-dr6-lens-values.ini b/examples/act-dr6-lens-values.ini index 4af5382a..10497b11 100644 --- a/examples/act-dr6-lens-values.ini +++ b/examples/act-dr6-lens-values.ini @@ -1,7 +1,7 @@ [cosmological_parameters] omch2 = 0.005 0.11 0.99 -ombh2 = 0.01 0.03 0.02233 +ombh2 = 0.01 0.02233 0.03 log1e10As = 1.61 3.0448 4.0 n_s = 0.86 0.96 1.06 cosmomc_theta = 0.005 0.0104 0.1 From 361becb9b3e9da976c8b419c9c9a687d7d0afdbf Mon Sep 17 00:00:00 2001 From: Joe Zuntz Date: Thu, 11 May 2023 17:00:26 +0100 Subject: [PATCH 04/12] update priors --- examples/act-dr6-lens-values.ini | 8 ++++---- examples/act-dr6-lens.ini | 15 ++++++++++----- 2 files changed, 14 insertions(+), 9 deletions(-) diff --git a/examples/act-dr6-lens-values.ini b/examples/act-dr6-lens-values.ini index 10497b11..820efbeb 100644 --- a/examples/act-dr6-lens-values.ini +++ b/examples/act-dr6-lens-values.ini @@ -1,17 +1,17 @@ [cosmological_parameters] -omch2 = 0.005 0.11 0.99 +omch2 = 0.08 0.12 0.16 ombh2 = 0.01 0.02233 0.03 -log1e10As = 1.61 3.0448 4.0 +log1e10As = 2.9 3.0448 3.2 n_s = 0.86 0.96 1.06 -cosmomc_theta = 0.005 0.0104 0.1 +cosmomc_theta = 0.009 0.0104 0.012 omega_k = 0.0 ;spatial curvature ;neutrinos mnu = 0.06 nnu = 3.046 -num_massive_neutrinos = 1 +num_massive_neutrinos = 3 ;helium yhe = 0.245341 ;helium mass fraction diff --git a/examples/act-dr6-lens.ini b/examples/act-dr6-lens.ini index cc7135b1..6beae73a 100644 --- a/examples/act-dr6-lens.ini +++ b/examples/act-dr6-lens.ini @@ -6,6 +6,10 @@ root = ${PWD} save_dir=output/act_dr6_lens fatal_errors=T +[nautilus] +n_live = 1500 +verbose = T + [pipeline] ; these names refer to sections later in the file: @@ -13,14 +17,15 @@ modules = consistency camb act_dr6_lens values = examples/act-dr6-lens-values.ini priors = examples/act-dr6-lens-priors.ini quiet=F -debug=T +debug=F timing=F [act_dr6_lens] file = ./likelihood/act-dr6-lens/act_dr6_lenslike_interface.py - +[output] +filename = output/act-dr6-lens.txt ; The consistency module translates between our chosen parameterization ; and any other that modules in the pipeline may want (e.g. camb) @@ -34,12 +39,12 @@ mode = cmb lmax = 4000 ;max ell to use for cmb calculation lens_margin = 1250 lens_potential_accuracy = 4 -feedback=2 ;amount of output to print +feedback=0 ;amount of output to print AccuracyBoost=1.0 ;CAMB accuracy boost parameter lSampleBoost = 1.0 lAccuracyBoost = 1.0 -do_tensors = True ;include tensor modes -do_lensing = true ;lensing is required w/ Planck data +do_tensors = T ;include tensor modes +do_lensing = T ;lensing is required w/ Planck data NonLinear = lens theta_H0_range = "20 100" halofit_version = mead2016 From 7a35e5f2f42ce293650ecb577f8d8548bbca2115 Mon Sep 17 00:00:00 2001 From: Joe Zuntz Date: Tue, 23 May 2023 12:22:12 +0100 Subject: [PATCH 05/12] add resume --- examples/act-dr6-lens.ini | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/act-dr6-lens.ini b/examples/act-dr6-lens.ini index 6beae73a..f6c2f1fb 100644 --- a/examples/act-dr6-lens.ini +++ b/examples/act-dr6-lens.ini @@ -1,6 +1,7 @@ [runtime] sampler = test root = ${PWD} +resume = T [test] save_dir=output/act_dr6_lens From 1b46d4d299ed031f079b0050f223bfecef2fc95a Mon Sep 17 00:00:00 2001 From: Joe Zuntz Date: Tue, 23 May 2023 12:23:44 +0100 Subject: [PATCH 06/12] comment update --- likelihood/act-dr6-lens/act_dr6_lenslike_interface.py | 1 + 1 file changed, 1 insertion(+) diff --git a/likelihood/act-dr6-lens/act_dr6_lenslike_interface.py b/likelihood/act-dr6-lens/act_dr6_lenslike_interface.py index 1349e9a9..3a87acc8 100644 --- a/likelihood/act-dr6-lens/act_dr6_lenslike_interface.py +++ b/likelihood/act-dr6-lens/act_dr6_lenslike_interface.py @@ -30,6 +30,7 @@ def execute(block, config): cl_te = block[names.cmb_cl, 'te'] / f1 cl_bb = block[names.cmb_cl, 'bb'] / f1 + # Slightly different normalization here f2 = ell * (ell + 1) cl_pp = block[names.cmb_cl, 'pp'] / f2 From 9d837c3c8339e123cd6d3d4ce5f31bfb48b6ddf4 Mon Sep 17 00:00:00 2001 From: Joe Zuntz Date: Thu, 13 Jul 2023 09:43:45 +0100 Subject: [PATCH 07/12] trying more options in act lensing test --- examples/act-dr6-lens.ini | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/examples/act-dr6-lens.ini b/examples/act-dr6-lens.ini index f6c2f1fb..29db48c9 100644 --- a/examples/act-dr6-lens.ini +++ b/examples/act-dr6-lens.ini @@ -11,6 +11,11 @@ fatal_errors=T n_live = 1500 verbose = T +[emcee] +walkers = 64 +samples = 1000 +nsteps = 5 + [pipeline] ; these names refer to sections later in the file: @@ -26,7 +31,7 @@ timing=F file = ./likelihood/act-dr6-lens/act_dr6_lenslike_interface.py [output] -filename = output/act-dr6-lens.txt +filename = output/act-dr6-lens-emcee.txt ; The consistency module translates between our chosen parameterization ; and any other that modules in the pipeline may want (e.g. camb) @@ -48,4 +53,4 @@ do_tensors = T ;include tensor modes do_lensing = T ;lensing is required w/ Planck data NonLinear = lens theta_H0_range = "20 100" -halofit_version = mead2016 +halofit_version = takahashi From b5b8a58b6cb893a93738d7aba773cbcb2c110e9e Mon Sep 17 00:00:00 2001 From: Joe Zuntz Date: Fri, 14 Jul 2023 22:14:38 +0100 Subject: [PATCH 08/12] Fix ACT PP normalization, thanks to David Dzingeleski --- examples/act-dr6-lens.ini | 2 ++ likelihood/act-dr6-lens/act_dr6_lenslike.py | 2 +- .../act_dr6_lenslike_interface.py | 19 ++++++++++++++----- 3 files changed, 17 insertions(+), 6 deletions(-) diff --git a/examples/act-dr6-lens.ini b/examples/act-dr6-lens.ini index 6beae73a..75c0a446 100644 --- a/examples/act-dr6-lens.ini +++ b/examples/act-dr6-lens.ini @@ -23,6 +23,8 @@ timing=F [act_dr6_lens] file = ./likelihood/act-dr6-lens/act_dr6_lenslike_interface.py +like_corrections = T +variant = actplanck_baseline [output] filename = output/act-dr6-lens.txt diff --git a/likelihood/act-dr6-lens/act_dr6_lenslike.py b/likelihood/act-dr6-lens/act_dr6_lenslike.py index 222c69e0..077408c0 100644 --- a/likelihood/act-dr6-lens/act_dr6_lenslike.py +++ b/likelihood/act-dr6-lens/act_dr6_lenslike.py @@ -345,7 +345,7 @@ def generic_lnlike(data_dict,ell_kk,cl_kk,ell_cmb,cl_tt,cl_ee,cl_te,cl_bb,trim_l bclkk = np.append(bclkk, d['binmat_planck'] @ clkk_planck) delta = d['data_binned_clkk'] - bclkk lnlike = -0.5 * np.dot(delta,np.dot(cinv,delta)) - return lnlike + return lnlike, bclkk # ================= diff --git a/likelihood/act-dr6-lens/act_dr6_lenslike_interface.py b/likelihood/act-dr6-lens/act_dr6_lenslike_interface.py index 3a87acc8..c53b0928 100644 --- a/likelihood/act-dr6-lens/act_dr6_lenslike_interface.py +++ b/likelihood/act-dr6-lens/act_dr6_lenslike_interface.py @@ -6,6 +6,7 @@ def setup(options): variant = options.get_string(option_section, 'variant', default='act_baseline') lens_only = options.get_bool(option_section, 'lens_only', default=False) like_corrections = options.get_bool(option_section, 'like_corrections', default=False) + like_only = options.get_bool(option_section, 'like_only', default=False) # lens_only use True if not combining with any primary CMB data # like_corrections should be False if lens_only is True @@ -13,6 +14,8 @@ def setup(options): # (covariance matrix) and `binmat_act` (binning matrix to be applied to a theory # curve starting at ell=0). data_dict = act_dr6_lenslike.load_data(variant,lens_only=lens_only,like_corrections=like_corrections) + + data_dict['cosmosis_like_only'] = like_only return data_dict @@ -24,20 +27,26 @@ def execute(block, config): # as well as the TT, EE, TE, BB CMB spectra (needed for likelihood corrections) # in uK^2 units. All of these are C_ell (not D_ell), no ell or 2pi factors. ell = block[names.cmb_cl, 'ell'] - f1 = ell * (ell + 1) / 2 / np.pi + f1 = ell * (ell + 1) / (2 * np.pi) cl_tt = block[names.cmb_cl, 'tt'] / f1 cl_ee = block[names.cmb_cl, 'ee'] / f1 cl_te = block[names.cmb_cl, 'te'] / f1 cl_bb = block[names.cmb_cl, 'bb'] / f1 # Slightly different normalization here - f2 = ell * (ell + 1) - cl_pp = block[names.cmb_cl, 'pp'] / f2 + cl_pp = block[names.cmb_cl, 'pp'] / f1 cl_kk = act_dr6_lenslike.pp_to_kk(cl_pp, ell) - # Then call - lnlike = act_dr6_lenslike.generic_lnlike(data_dict,ell, cl_kk, ell, cl_tt, cl_ee, cl_te, cl_bb) + # Then call the act code + lnlike, bclkk = act_dr6_lenslike.generic_lnlike(data_dict,ell, cl_kk, ell, cl_tt, cl_ee, cl_te, cl_bb) block[names.likelihoods, 'act_dr6_lens_like'] = lnlike + if not data_dict['cosmosis_like_only']: + block[names.data_vector, 'act_dr6_lens_theory'] = bclkk + block[names.data_vector, 'act_dr6_lens_data'] = data_dict['data_binned_clkk'] + block[names.data_vector, 'act_dr6_lens_covariance'] = data_dict['cov'] + block[names.data_vector, 'act_dr6_lens_inverse_covariance'] = data_dict['cinv'] + + return 0 \ No newline at end of file From 0c4c6323e162c51b76f59125d409d66aa01042e4 Mon Sep 17 00:00:00 2001 From: Joe Zuntz Date: Mon, 17 Jul 2023 09:35:02 +0100 Subject: [PATCH 09/12] save extra output in act example --- examples/act-dr6-lens.ini | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/act-dr6-lens.ini b/examples/act-dr6-lens.ini index 53ac5b16..f47b97ab 100644 --- a/examples/act-dr6-lens.ini +++ b/examples/act-dr6-lens.ini @@ -25,7 +25,7 @@ priors = examples/act-dr6-lens-priors.ini quiet=F debug=F timing=F - +extra_output=cosmological_parameters/sigma_8 cosmological_parameters/omega_m [act_dr6_lens] file = ./likelihood/act-dr6-lens/act_dr6_lenslike_interface.py @@ -33,7 +33,7 @@ like_corrections = T variant = actplanck_baseline [output] -filename = output/act-dr6-lens-emcee.txt +filename = output/act-dr6-lens.txt ; The consistency module translates between our chosen parameterization ; and any other that modules in the pipeline may want (e.g. camb) From 6613e0b2fb46b74b4ab17ca5709d9f244d304973 Mon Sep 17 00:00:00 2001 From: Joe Zuntz Date: Mon, 17 Jul 2023 10:39:06 +0100 Subject: [PATCH 10/12] Add ACT DR6 to CI --- .github/workflows/ci.yml | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index fa015a4f..8d91e817 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -169,6 +169,13 @@ jobs: cosmosis examples/des-y3.ini -v halo_model_parameters.logT_AGN=8.2 -p camb.halofit_version=mead2020_feedback | tee output/des-y3-mead.log grep -e 'Likelihood = 6049.94' -e 'Likelihood = 6049.00' output/des-y3-mead.log + - name: ACT DR6 Lensing + shell: bash -l {0} + run: | + source cosmosis-configure + cosmosis examples/act-dr6-lens.ini | tee output/act-dr6.log + grep -e "Likelihood = -9.89" -e Likelihood = -9.90" output/act-dr6.log + - name: DES Y3 6x2pt shell: bash -l {0} run: | From e4a33f3dca610f7b2c8ece4e04805846666d8ce7 Mon Sep 17 00:00:00 2001 From: Joe Zuntz Date: Mon, 17 Jul 2023 10:51:47 +0100 Subject: [PATCH 11/12] add data download and cache --- .github/workflows/ci.yml | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 8d91e817..7c642768 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -169,6 +169,19 @@ jobs: cosmosis examples/des-y3.ini -v halo_model_parameters.logT_AGN=8.2 -p camb.halofit_version=mead2020_feedback | tee output/des-y3-mead.log grep -e 'Likelihood = 6049.94' -e 'Likelihood = 6049.00' output/des-y3-mead.log + - uses: actions/cache@v2 + name: ACT Data Cache + id: cache-act + with: + path: likelihood/act-dr6-lens/data/v1.1 + key: ${{ runner.os }}-act-dr6-v1.1 + + - name: Download ACT DR6 Lensing Data + if: steps.cache-act.outputs.cache-hit != 'true' + run: | + cd likelihood/act-dr6-lens/ + ./get-act-data.sh + - name: ACT DR6 Lensing shell: bash -l {0} run: | From 13b596ed1c3982217c15592148fca17ab7e8ad6e Mon Sep 17 00:00:00 2001 From: Joe Zuntz Date: Mon, 17 Jul 2023 11:09:19 +0100 Subject: [PATCH 12/12] fix CI syntax --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 7c642768..e14eb363 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -187,7 +187,7 @@ jobs: run: | source cosmosis-configure cosmosis examples/act-dr6-lens.ini | tee output/act-dr6.log - grep -e "Likelihood = -9.89" -e Likelihood = -9.90" output/act-dr6.log + grep -e 'Likelihood = -9.89' -e 'Likelihood = -9.90' output/act-dr6.log - name: DES Y3 6x2pt shell: bash -l {0}