From 44e46c8b78412821f8160f89dd7e0a77f4224773 Mon Sep 17 00:00:00 2001 From: Eunseong Lee <36232273+eunseongleee@users.noreply.github.com> Date: Wed, 13 Aug 2025 15:46:51 -0400 Subject: [PATCH 01/14] frequency null test updated using an ILC filter that's used for the baseline; comparison with public DR6 lensing is added; tiny update for reading both SZ and optical sample catalogue --- stack.py | 170 +++++++++++++++++++++++++++++++++++++++---------------- utils.py | 124 +++++++++++++++++----------------------- 2 files changed, 175 insertions(+), 119 deletions(-) diff --git a/stack.py b/stack.py index c404279..a60c781 100644 --- a/stack.py +++ b/stack.py @@ -7,8 +7,8 @@ matplotlib.use('Agg') import matplotlib.pyplot as plt import utils as cutils -from pixell import enmap, reproject, enplot, utils, wcsutils -from orphics import maps, mpi, io, stats, cosmology +from pixell import enmap, reproject, utils, wcsutils, curvedsky +from orphics import maps, mpi, io, stats, cosmology, lensing from scipy.optimize import curve_fit from numpy import save import time @@ -90,11 +90,6 @@ # reproject the Planck map (healpix -> CAR) fshape, fwcs = enmap.fullsky_geometry(res=2.0 * utils.arcmin, proj="car") - # this doesn't work with the latest pixell version (v0.20.3) - # pmap = reproject.enmap_from_healpix( - # plc_map, fshape, fwcs, ncomp=1, unit=1, lmax=6000, rot="gal,equ" - # ) - # reading the input map p_map = np.atleast_2d(hp.read_map(plc_map, field=tuple(range(0,1)))).astype(np.float64) @@ -115,7 +110,6 @@ else: act_map = (paths.act_data + data_choice.hres) famap_150 = enmap.read_map(act_map, delayed=False) - print(np.shape(famap_150)) # SZ cluster model image subtraction for 150 GHz if args.hres_grad: @@ -128,7 +122,6 @@ amap_150 = famap_150 - enmap.read_map(paths.act_data + data_choice.hres_model_150) else: amap_150 = famap_150 - # ACT 90 GHz coadd map if not(args.full_sim_index is None): @@ -140,7 +133,7 @@ # SZ cluster model image subtraction for 90 GHz if args.hres_grad: if not (args.grad_noszsub): - gamap_090 = famap_90 - enmap.read_map(paths.act_data + data_choice.hres_model_090) + gamap_90 = famap_90 - enmap.read_map(paths.act_data + data_choice.hres_model_090) else: gamap_90 = famap_90 @@ -156,11 +149,11 @@ assert not(args.no_150) assert not(args.rand_rot) act_map = ( - paths.coadd_data + f"{tags.apstr}_s08_{tags.s19str}_cmb_f150_night_srcfree_map.fits" + paths.act_data + data_choice.hres_150_night ) namap_150 = enmap.read_map(act_map, delayed=False, sel=np.s_[0, ...]) act_map = ( - paths.coadd_data + f"{tags.apstr}_s08_{tags.s19str}_cmb_f090_night_srcfree_map.fits" + paths.act_data + data_choice.hres_090_night ) namap_90 = enmap.read_map(act_map, delayed=False, sel=np.s_[0, ...]) @@ -170,18 +163,25 @@ # Inv var map for 90 GHz ivar_map = (paths.act_data + data_choice.hres_ivar) - # if data_choice.hres_map == 'dr6': imap_90 = enmap.read_map(ivar_map, delayed=False) - # else: imap_90 = enmap.read_map(ivar_map, delayed=False, sel=np.s_[0, ...]) ##### fix here - try: imap_90 = enmap.read_map(ivar_map, delayed=False, sel=np.s_[0, ...]) except: imap_90 = enmap.read_map(ivar_map, delayed=False) + rms_map = maps.rms_from_ivar( imap_90, cylindrical=True ) # convert to RMS noise map + + if args.dr6_lensing: + kappa_map = paths.act_data + "release/dr6_lensing_v1/maps/baseline/kappa_alm_data_act_dr6_lensing_v1_baseline.fits" + kmap = np.nan_to_num(hp.read_alm(kappa_map).astype(np.complex128)) + # reading mask for geometry, mask is NOT applied to the map + mask = enmap.read_map(paths.act_data + "DR6_lensing/masks/act_mask_20220316_GAL060_rms_70.00_d2sk.fits") + k_map = curvedsky.alm2map(kmap, enmap.empty(mask.shape, mask.wcs, dtype=np.float64)) + + # stamp size and resolution stamp_width_deg = args.swidth / 60.0 # stamp_width_arcmin: 128.0 pixel = args.pwidth # pix_width_arcmin: 0.5 @@ -377,6 +377,13 @@ def fit_p1d( bin_edges = np.arange(0, args.arcmax, args.arcstep) # 15 arcmin, 1.5 arcmin centers = (bin_edges[1:] + bin_edges[:-1]) / 2.0 +# for kappa L +ell_edges = np.arange(klmin, klmax, 200) +ell_cents = (ell_edges[1:] + ell_edges[:-1]) / 2.0 + + + + """ !! BINNING """ @@ -415,6 +422,25 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): return ret, tret + + +#------------------------------------------------------------------------------- + +# agora sample mean +M200c = 6.37e14 +z = 0.6 + +rstamp = 128.8 * utils.arcmin # args.swidth = 128 +res = pixel * utils.arcmin +rmin = 0. +rmax = 10 * utils.arcmin +rwidth = 1.0 * utils.arcmin + +_,_,_,_,_,_,_,_,k2dmap,_ = lensing.kappa_nfw_profiley(mass=M200c,conc=None,z=z,z_s=1100.,background='critical',delta=200,apply_filter=True,lmin=klmin,lmax=klmax,res=res,rstamp=rstamp,rmin=rmin,rmax=rmax,rwidth=rwidth) + +#------------------------------------------------------------------------------- + + """ !! LOOP OVER ASSIGNED TASKS """ @@ -619,6 +645,20 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): else: pstamp, astamp_150, astamp_90 = csim.get_obs(task) + if args.dr6_lensing: + # cut out a stamp from the ACT DR6 reconstructed lensing map + kstamp = reproject.thumbnails( + k_map, + coords, + r=maxr, + res=pixel * utils.arcmin, + proj="tan", + oversample=2, + pixwin=False + ) + + + """ !! COSINE TAPER """ @@ -637,9 +677,6 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): act_stamp_150 = astamp_150 * taper act_stamp_90 = astamp_90 * taper - if args.freq_null: - act_stamp_fnull = (astamp_150 - astamp_90) * taper - if not (args.hres_grad): plc_stamp = pstamp * taper else: @@ -700,6 +737,8 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): if args.no_filter: s.add_to_stack('a150_cmb', astamp_150*sweight) s.add_to_stack('a90_cmb', astamp_90*sweight) + # s.add_to_stack('a150_cmb', astamp_150) + # s.add_to_stack('a90_cmb', astamp_90) sz150 = bin(astamp_150, modrmap * (180 * 60 / np.pi), bin_edges) sz150w = bin(astamp_150*sweight, modrmap * (180 * 60 / np.pi), bin_edges) sz90 = bin(astamp_90, modrmap * (180 * 60 / np.pi), bin_edges) @@ -777,6 +816,7 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): # j = j + 1 # continue # commented for now + """ !! STAMP FFTs """ @@ -787,9 +827,6 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): if args.day_null: nk150 = enmap.fft(nact_stamp_150, normalize="phys") nk90 = enmap.fft(nact_stamp_90, normalize="phys") - - if args.freq_null: - fnk = enmap.fft(act_stamp_fnull, normalize="phys") if not (args.hres_grad): kp = enmap.fft(plc_stamp, normalize="phys") @@ -863,9 +900,6 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): act_cents, act_p1d_90 = lbinner.bin(pow(k90, k90) / w2) act_cents, act_p1d_150_90 = lbinner.bin(pow(k150, k90) / w2) - if args.freq_null: - act_cents, act_p1d = lbinner.bin(pow(fnk, fnk) / w2) - if not (args.hres_grad): plc_cents, plc_p1d = lbinner.bin(pow(kp, kp) / w2) else: @@ -912,21 +946,7 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): rms=0, lmin=defaults.highres_fit_ellmin, lmax=defaults.highres_fit_ellmax, - ) - - if args.freq_null: - tclaa_fn = fit_p1d( - l_edges, - act_cents, - act_p1d, - "act", - modlmap, - bfunc150, - bfunc150, - rms=defaults.highres_fiducial_rms, - lmin=defaults.highres_fit_ellmin, - lmax=defaults.highres_fit_ellmax, - ) + ) if not (args.hres_grad): tclpp = fit_p1d( @@ -1024,7 +1044,6 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): act_kmap150 = k150 / act_150_kbeam2d act_kmap90 = k90 / act_90_kbeam2d act_kmap = act_kmap150 - act_kmap90 - tclaa = tclaa_fn / (act_150_kbeam2d ** 2.0) if args.no_150: act_kmap = k90 / act_90_kbeam2d @@ -1116,6 +1135,49 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): # transform to real space for unweighted stack kappa = enmap.ifft(krecon, normalize="phys").real + + + + #--------------------------------------------------------------------------- + ktemp = k2dmap # this is tapered and filtered template (in real space) + + template = enmap.fft(ktemp, normalize="phys") + assert np.all(np.isfinite(template)) + + p2d_cross_r = (krecon * template.conj()).real + p2d_auto_template = (template * template.conj()).real + + p1d_cross_r = bin(p2d_cross_r, modlmap, ell_edges) + p1d_auto_template = bin(p2d_auto_template, modlmap, ell_edges) + + s.add_to_stats("p1d_cross_r", p1d_cross_r) + s.add_to_stats("p1d_auto_template", p1d_auto_template) + + s.add_to_stack("temp", k2dmap) + binned_temp = bin(k2dmap, modrmap * (180 * 60 / np.pi), bin_edges) + s.add_to_stats("temp1d", binned_temp) + #--------------------------------------------------------------------------- + + + #--------------------------------------------------------------------------- + # for DR6 lensing stack + if args.dr6_lensing: + + k_stamp = maps.filter_map(kstamp * taper, kmask) # this is tapered and filtered DR6 kappa (real space) + + ikmap = enmap.fft(k_stamp, normalize="phys") + assert np.all(np.isfinite(ikmap)) + + p2d_cross_k = (ikmap * template.conj()).real # cross with template + p1d_cross_k = bin(p2d_cross_k, modlmap, ell_edges) + s.add_to_stats("p1d_cross_k", p1d_cross_k) + + s.add_to_stack("kstamp", k_stamp) + binned_kappa = bin(k_stamp, modrmap * (180 * 60 / np.pi), bin_edges) + s.add_to_stats("tk1d", binned_kappa) + #--------------------------------------------------------------------------- + + """ !! REJECT WEIRD KAPPA """ @@ -1228,9 +1290,10 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): if not(args.is_meanfield) and not(args.inject_sim): s.add_to_stats("data", (z,weight,*[cdata[key][i] for key in sorted(cdata.keys())])) s.add_to_stats("redshift", (z,)) - #s.add_to_stats("mass", (cdata['mass'][i],)) - #s.add_to_stats("y0s", (cdata['y0s'][i],)) - #s.add_to_stats("wmass", (cdata['mass'][i] * weight,)) + try: + s.add_to_stats("mass", (cdata['mass'][i],)) + except: + pass j = j + 1 @@ -1356,7 +1419,11 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): np.savetxt(f"{paths.savedir}/bin_edges.txt", bin_edges) if not(args.is_meanfield) and not (args.debug_stack): np.savetxt(f"{paths.savedir}/profiles.txt",s.vectors['k1d']) - #np.savetxt(f"{paths.savedir}/z_mass_y.txt", np.c_[s.vectors['redshift'], s.vectors['mass'], s.vectors['y0s']]) + try: + np.savetxt(f"{paths.savedir}/z_mass.txt", np.c_[s.vectors['redshift'], s.vectors['mass']]) + # np.savetxt(f"{paths.savedir}/z_mass_y.txt", np.c_[s.vectors['redshift'], s.vectors['mass'], s.vectors['y0s']]) + except: + pass for ctkey in [ "selected" @@ -1373,10 +1440,17 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): assert N_stamp == s.stack_count["kmap"] assert N_stamp == s.vectors["kw"].shape[0] - # if not (args.is_meanfield): - # kappa_w = s.vectors['kw'].sum() - # w_mass = s.vectors['wmass'].sum(axis=0)/kappa_w - # print("weighted SZ mean mass : ", w_mass, "1e14 Msun") + if not (args.is_meanfield) and not (args.debug_stack): + print("mean redshift:", s.vectors['redshift'].sum() / N_stamp) + try: + print("mean SZ mass [M500c]:", s.vectors['mass'].sum() / N_stamp) + except: + pass + + # if not (args.is_meanfield): + # kappa_w = s.vectors['kw'].sum() + # w_mass = s.vectors['wmass'].sum(axis=0)/kappa_w + # print("weighted SZ mean mass : ", w_mass, "1e14 Msun") # Some debug plots if requested if args.debug_powers: diff --git a/utils.py b/utils.py index fecaada..fddb13c 100644 --- a/utils.py +++ b/utils.py @@ -238,6 +238,11 @@ def initialize_pipeline_config(): parser.add_argument( "--decmin", type=float, default=None, help="Minimum declination in degree." ) + parser.add_argument( + "--dr6-lensing", + action="store_true", + help="Stack dr6 reconstructed lensing (no reconstruction only stacking!)." + ) args = parser.parse_args() if args.hres_lmin is None: @@ -312,57 +317,40 @@ def initialize_pipeline_config(): paths.savedir = savedir return start_time,paths,defaults,args,tags,rank,data_choice -def cut_z_sn(ras,decs,sns,zs,zmin,zmax,snmin,snmax,y0s,y0min,y0max,decmin,mass): +def cut_z_sn(ras,decs,sns,zs,zmin=None,zmax=None,snmin=None,snmax=None,y0s=None,y0min=None,y0max=None,decmin=None,mass=None): + + def select_only(mask): + nonlocal ras, decs, sns, zs, y0s, mass + ras = ras[mask] + decs = decs[mask] + sns = sns[mask] + zs = zs[mask] + if y0s is not None: + y0s = y0s[mask] + if mass is not None: + mass = mass[mask] + if zmin is not None: - ras = ras[zs>zmin] - decs = decs[zs>zmin] - sns = sns[zs>zmin] - y0s = y0s[zs>zmin] - mass = mass[zs>zmin] - zs = zs[zs>zmin] + select_only(zs > zmin) if zmax is not None: - ras = ras[zs<=zmax] - decs = decs[zs<=zmax] - sns = sns[zs<=zmax] - y0s = y0s[zs<=zmax] - mass = mass[zs<=zmax] - zs = zs[zs<=zmax] + select_only(zs <= zmax) + if snmin is not None: - ras = ras[sns>snmin] - decs = decs[sns>snmin] - zs = zs[sns>snmin] - y0s = y0s[sns>snmin] - mass = mass[sns>snmin] - sns = sns[sns>snmin] + select_only(sns > snmin) if snmax is not None: - ras = ras[sns<=snmax] - decs = decs[sns<=snmax] - zs = zs[sns<=snmax] - y0s = y0s[sns<=snmax] - mass = mass[sns<=snmax] - sns = sns[sns<=snmax] - if y0min is not None: - ras = ras[y0s>y0min] - decs = decs[y0s>y0min] - zs = zs[y0s>y0min] - sns = sns[y0s>y0min] - mass = mass[y0s>y0min] - y0s = y0s[y0s>y0min] - if y0max is not None: - ras = ras[y0s<=y0max] - decs = decs[y0s<=y0max] - zs = zs[y0s<=y0max] - sns = sns[y0s<=y0max] - mass = mass[y0s<=y0max] - y0s = y0s[y0s<=y0max] + select_only(sns <= snmax) + + if y0s is not None: + if y0min is not None: + select_only(y0s > y0min) + if y0max is not None: + select_only(y0s <= y0max) + if decmin is not None: - ras = ras[np.abs(decs)<=decmin] - zs = zs[np.abs(decs)<=decmin] - y0s = y0s[np.abs(decs)<=decmin] - sns = sns[np.abs(decs)<=decmin] - mass = mass[np.abs(decs)<=decmin] - decs = decs[np.abs(decs)<=decmin] - return ras,decs,sns,zs,y0s,mass + select_only(np.abs(decs) <= decmin) + + return ras, decs, sns, zs, y0s, mass + def catalog_interface(cat_type,is_meanfield,nmax=None,zmin=None,zmax=None,bcg=False,snmin=None,snmax=None,y0min=None,y0max=None,decmin=None): data = {} @@ -384,6 +372,7 @@ def catalog_interface(cat_type,is_meanfield,nmax=None,zmin=None,zmax=None,bcg=Fa sns = hdu[1].data['SNR'][ras>=0] y0s = hdu[1].data['fixed_y_c'][ras>=0] mass = hdu[1].data['M500'][ras>=0] + # mass = hdu[1].data['M200m'][ras>=0] ras = ras[ras>=0] else: ras = hdu[1].data['RADeg'] @@ -555,16 +544,16 @@ def catalog_interface(cat_type,is_meanfield,nmax=None,zmin=None,zmax=None,bcg=Fa elif cat_type=='des_redmapper': if is_meanfield: - catalogue_name = paths.data+ 'y3_gold_2.2.1_wide_sofcol_run_redmapper_v6.4.22_randcat_z0.10-0.95_lgt020_vl02.fit' + catalogue_name = paths.des_data+ 'y3_gold_2.2.1_wide_sofcol_run2_redmapper_v6.4.22+2_randcat_z0.10-0.95_lgt020_vl02.fit' else: - catalogue_name = paths.data+ 'y3_gold_2.2.1_wide_sofcol_run_redmapper_v6.4.22_lgt20_vl02_catalog.fit' + catalogue_name = paths.des_data+ 'y3_gold_2.2.1_wide_sofcol_run2_redmapper_v6.4.22+2_lgt20_vl02_catalog.fit' hdu = fits.open(catalogue_name) ras = hdu[1].data['RA'] decs = hdu[1].data['DEC'] zs = hdu[1].data['Z_LAMBDA' if not(is_meanfield) else 'ZTRUE'] sns = hdu[1].data['LAMBDA_CHISQ' if not(is_meanfield) else 'LAMBDA_IN'] - ras,decs,sns,zs = cut_z_sn(ras,decs,sns,zs,zmin,zmax,snmin,snmax) + ras,decs,sns,zs,_,_ = cut_z_sn(ras,decs,sns,zs,zmin,zmax,snmin,snmax) ras = ras[:nmax] decs = decs[:nmax] @@ -639,8 +628,7 @@ def catalog_interface(cat_type,is_meanfield,nmax=None,zmin=None,zmax=None,bcg=Fa zs = zs[sns>0] sns = sns[sns>0] - ras,decs,sns,zs = cut_z_sn(ras,decs,sns,zs,zmin,zmax,snmin,snmax) - + ras,decs,sns,zs,_,_ = cut_z_sn(ras,decs,sns,zs,zmin,zmax,snmin,snmax) sns = sns[:nmax] data['lams'] = sns @@ -659,7 +647,7 @@ def catalog_interface(cat_type,is_meanfield,nmax=None,zmin=None,zmax=None,bcg=Fa else: catalogue_name = paths.data+ 'camira_s19a_wide_sm_v1_01z11.dat' ras,decs,sns,zs = np.loadtxt(catalogue_name,unpack=True) - ras,decs,sns,zs = cut_z_sn(ras,decs,sns,zs,zmin,zmax,snmin,snmax) + ras,decs,sns,zs,_,_ = cut_z_sn(ras,decs,sns,zs,zmin,zmax,snmin,snmax) ras = ras[:nmax] decs = decs[:nmax] @@ -683,17 +671,14 @@ def catalog_interface(cat_type,is_meanfield,nmax=None,zmin=None,zmax=None,bcg=Fa data['lams'] = ws*0 - else: raise NotImplementedError return ras,decs,zs,ws,data def load_beam(freq): - #if freq=='f150': fname = paths.data+'s16_pa2_f150_nohwp_night_beam_tform_jitter.txt' - #elif freq=='f090': fname = paths.data+'s16_pa3_f090_nohwp_night_beam_tform_jitter.txt' - if freq=='f150': fname = paths.data+'corrected_beam_150.txt' - elif freq=='f090': fname = paths.data+'corrected_beam_090.txt' + if freq=='f150': fname = paths.act_data+'coadds/20240323_simple/beam_f150_tform.txt' + elif freq=='f090': fname = paths.act_data+'coadds/20240323_simple/beam_f090_tform.txt' ls,bls = np.loadtxt(fname,usecols=[0,1],unpack=True) assert ls[0]==0 bls = bls / bls[0] @@ -959,14 +944,14 @@ def load_vrec_catalog_boss(pathOutCatalog): def postprocess(stack_path,mf_path,save_name=None,ignore_param=False,args=None,ignore_last=None): - if mf_path!="": + if mf_path != "": smf_path = mf_path if (ignore_last is None) else mf_path[:-ignore_last] mf_paramstr = re.search(rf'plmin_(.*?)_meanfield', smf_path).group(1) sstack_path = stack_path if (ignore_last is None) else stack_path[:-ignore_last] st_paramstr = re.search(rf'plmin_(.*)', sstack_path).group(1) if not(ignore_param): - if mf_path!="": + if mf_path != "": try: assert mf_paramstr==st_paramstr except: @@ -989,10 +974,10 @@ def postprocess(stack_path,mf_path,save_name=None,ignore_param=False,args=None,i if data is not None: io.save_cols(f'{save_dir}/{save_name}_catalog_data.txt',[data[key] for key in sorted(data.keys())],header=' '.join([key for key in sorted(data.keys())])) - if mf_path!="": + if mf_path != "": s_mf, shape_mf, wcs_mf = load_dumped_stats(mf_path) - if mf_path!="": + if mf_path != "": assert np.all(shape_stack==shape_mf) assert wcsutils.equal(wcs_stack,wcs_mf) assert np.all(shape_stack==kmask.shape) @@ -1004,7 +989,7 @@ def postprocess(stack_path,mf_path,save_name=None,ignore_param=False,args=None,i crop = int(args.cwidth / defaults.pix_width_arcmin) unweighted_stack,nmean_weighted_kappa_stack,opt_weighted_kappa_stack,opt_binned,opt_covm,opt_corr,opt_errs,binned,covm,corr,errs = analyze(s_stack,wcs) - if mf_path!="": + if mf_path != "": mf_unweighted_stack,mf_nmean_weighted_kappa_stack,mf_opt_weighted_kappa_stack,mf_opt_binned,mf_opt_covm,mf_opt_corr,mf_opt_errs,mf_binned,mf_covm,mf_corr,mf_errs = analyze(s_mf,wcs) # if profs is not None: @@ -1031,7 +1016,7 @@ def postprocess(stack_path,mf_path,save_name=None,ignore_param=False,args=None,i plot(f"{save_dir}/{save_name}_unweighted_nomfsub.png",unweighted_stack,tap_per,pad_per,crop=None,lim=args.plim) plot(f"{save_dir}/{save_name}_unweighted_nomfsub_zoom.png",unweighted_stack,tap_per,pad_per,crop=crop,lim=args.plim) - if mf_path!="": + if mf_path != "": # Opt weighted stamp = opt_weighted_kappa_stack - mf_opt_weighted_kappa_stack plot(f"{save_dir}/{save_name}_opt_weighted_mfsub.png",stamp,tap_per,pad_per,crop=None,lim=args.plim) @@ -1096,7 +1081,7 @@ def postprocess(stack_path,mf_path,save_name=None,ignore_param=False,args=None,i io.plot_img(corr,f'{save_dir}/{save_name}_corr.png') pl = io.Plotter(xyscale='linlin', xlabel='$\\theta$ [arcmin]', ylabel='$\\kappa$') - if mf_path!="": + if mf_path != "": pl.add_err(cents, opt_binned - mf_opt_binned, yerr=opt_errs,ls="-",label="Filtered kappa, mean-field subtracted (optimal)") pl.add_err(cents+0.2, binned - mf_binned, yerr=errs,ls="-",label="Filtered kappa, mean-field subtracted") pl.add_err(cents, mf_opt_binned, yerr=mf_opt_errs,label="Mean-field (optimal)",ls="-",alpha=0.5) @@ -1111,7 +1096,7 @@ def postprocess(stack_path,mf_path,save_name=None,ignore_param=False,args=None,i pl.done(f'{save_dir}/{save_name}_profile.png') pl = io.Plotter(xyscale='linlin', xlabel='$\\theta$ [arcmin]', ylabel='$\\kappa$') - if mf_path!="": + if mf_path != "": pl.add_err(cents, opt_binned - mf_opt_binned, yerr=opt_errs,ls="-",label="Filtered kappa, mean-field subtracted (optimal)") pl.add_err(cents, mf_opt_binned, yerr=mf_opt_errs,label="Mean-field (optimal)",ls="-",alpha=0.5) else: @@ -1144,7 +1129,7 @@ def postprocess(stack_path,mf_path,save_name=None,ignore_param=False,args=None,i io.save_cols(f'{save_dir}/{save_name}_profile.txt', (cents, ret_data)) io.save_cols(f'{save_dir}/{save_name}_profile_errs.txt', (cents, errs)) io.save_cols(f'{save_dir}/{save_name}_mf.txt', (cents, mf_binned)) - if mf_path!="": + if mf_path != "": io.save_cols(f'{save_dir}/{save_name}_mf_errs.txt', (cents, mf_opt_errs)) np.savetxt(f'{save_dir}/{save_name}_opt_covm.txt', ret_opt_cov) np.savetxt(f'{save_dir}/{save_name}_covm.txt', ret_cov) @@ -1175,6 +1160,7 @@ def postprocess(stack_path,mf_path,save_name=None,ignore_param=False,args=None,i cov = opt_covm[:nbins,:nbins] fbin_edges = bin_edges[:nbins+1] fcents = cents[:nbins] + lnlikes,like_fit,fit_mass,mass_err,fprofiles,fit_profile =lensing.fit_nfw_profile(profile,cov,masses,z,conc,cc,shape,wcs,fbin_edges,lmax=0,lmin=0,overdensity=args.overdensity,critical=args.critical,at_cluster_z=args.at_z0,mass_guess=mguess,sigma_guess=merr_guess,kmask=kmask,sigma_mis=sigma_mis) print("Fit mass : " , fit_mass/1e14,mass_err/1e14) @@ -1226,8 +1212,4 @@ def postprocess(stack_path,mf_path,save_name=None,ignore_param=False,args=None,i pl.hline(y=1) pl.done(f'{save_dir}/{save_name}_theory_comp_ratio.png') - - - return cents,ret_opt_data,ret_opt_cov - - + return cents,ret_opt_data,ret_opt_cov \ No newline at end of file From 59f9ec649b9150b46be02f867f0e4d3f5ef21a42 Mon Sep 17 00:00:00 2001 From: Eunseong Lee <36232273+eunseongleee@users.noreply.github.com> Date: Thu, 14 Aug 2025 14:39:26 -0400 Subject: [PATCH 02/14] attempt to separate the templates for lensing and cross-correlation --- analysis.py | 14 ++++++++++---- choices.yaml | 6 ++++-- 2 files changed, 14 insertions(+), 6 deletions(-) diff --git a/analysis.py b/analysis.py index ea20e22..6da5e6a 100644 --- a/analysis.py +++ b/analysis.py @@ -39,12 +39,18 @@ def __init__(self, self.debug = debug - # Load template + # Load template for simulation thetas,kappa_1h,kappa_2h = np.loadtxt(c.template_file,unpack=True) kappa = kappa_1h + kappa_2h self.thetas = thetas self.template_kappa_1d = kappa + # Load template for cross-correlation + thetas_cross,kappa_1h_cross,kappa_2h_cross = np.loadtxt(c.cross_template_file,unpack=True) + kappa_cross = kappa_1h_cross + kappa_2h_cross + self.thetas_cross = thetas_cross + self.template_kappa_1d_cross = kappa_cross + if atype=='flatsky_sim': # If we are doing a padded high-res flat sim # (e.g. for multiplicative bias correction) @@ -66,7 +72,7 @@ def __init__(self, self.rbinner = stats.bin2D(self.modrmap,self.rbin_edges) self.rcents = self.rbinner.cents self.modlmap = enmap.modlmap(self.shape,self.wcs) - self.template = enmap.enmap(maps.interp(thetas,kappa)(self.modrmap),self.wcs) + self.template_cross = enmap.enmap(maps.interp(thetas_cross,kappa_cross)(self.modrmap),self.wcs) self.reconstructor = Recon(self.shape,self.wcs, @@ -97,7 +103,7 @@ def __init__(self, self.Lmin = c.Lmin self.Lmax = c.Lmax - self.ktemplate = enmap.fft(self.template * self.ctaper,normalize='phys') + self.ktemplate = enmap.fft(self.template_cross * self.ctaper,normalize='phys') cents,p1d,p2d = self.power(self.ktemplate,self.ktemplate) self.template_auto_p1d = p1d if debug and self.rank==0: @@ -439,4 +445,4 @@ def recon(self,imapx,imapy=None,p2d_plot=None): # Fourier space lens reconstruction krecon = self.cqe.reconstruct(self.feed_dict, xname="X_l1", yname="Y_l2", physical_units=True) - return krecon + return krecon \ No newline at end of file diff --git a/choices.yaml b/choices.yaml index b1efdb0..e8b5ff9 100644 --- a/choices.yaml +++ b/choices.yaml @@ -56,7 +56,9 @@ Lmax: 5000 # This is the template used in all cross-correlations # You don't need to know how it was made as long as the data and fits are always done with the same one # This particular one is a 1h+2h NFW mass=2e14 at z=1 with some power law extrapolations -template_file: data/kappa_profile.txt +cross_template_file: data/kappa_profile.txt +# This is the template used to simulate the cluster +# This particular one is a 1h+2h NFW mass=1e14 at z=1 with some power law extrapolations +template_file: data/kappa_profile_1e14.txt # For real-space fits, maximum radius to which to go out to (which determines the size of the covmat) theta_max_real_arc: 10.0 - From ce569feadd9023714ee8650111365cf842e7a7ec Mon Sep 17 00:00:00 2001 From: Eunseong Lee <36232273+eunseongleee@users.noreply.github.com> Date: Fri, 24 Oct 2025 13:35:17 -0400 Subject: [PATCH 03/14] ILC of 150 and 90 GHz included --- recon_sim.py | 783 +++++++++++++++++++++++++++++++++++++++++---------- 1 file changed, 633 insertions(+), 150 deletions(-) diff --git a/recon_sim.py b/recon_sim.py index 2ac2bb1..8d14abf 100644 --- a/recon_sim.py +++ b/recon_sim.py @@ -2,12 +2,12 @@ from astropy.io import fits import matplotlib matplotlib.use("Agg") -import matplotlib.pyplot as plt import numpy as np from pixell import enmap, reproject, utils, wcsutils, bunch -from orphics import mpi, maps, stats, io, cosmology -from symlens import qe +from orphics import mpi, maps, stats, io, cosmology, lensing +import symlens import argparse +import sys, os start = t.time() @@ -28,13 +28,10 @@ "which_cat", type=str, help="Choose the catalogue type e.g. halo or tsz." ) parser.add_argument( - "hres_choice", type=str, help="Choose the map for high resolution leg e.g. cmb, cmb_tsz, etc." + "--cmb-tsz", action="store_true", help="CMB + tSZ map for high resolution leg." ) parser.add_argument( - "grad_choice", type=str, help="Choose the map for gradient leg e.g. cmb, cmb_tsz, etc." -) -parser.add_argument( - "--is-observed", action="store_true", help="Beam and noise have been applied to the map." + "--cmb-msub", action="store_true", help="CMB + tSZ and model image subtraction for high resolution leg." ) parser.add_argument( "--is-meanfield", action="store_true", help="This is a mean-field run." @@ -43,17 +40,26 @@ "--is-test", action="store_true", help="This is a test run for first 10 entries." ) parser.add_argument( - "--full-sample", action="store_true", help="Entire sample will be used (SNR > 4)." + "--full-sample", action="store_true", help="Entire sample with SNR > 4 will be used.(Default cut is SNR > 5.5)" +) +parser.add_argument( + "--highsnr-sample", action="store_true", help="Higher SNR cut sample (SNR > 7) will be used." +) +parser.add_argument( + "--snmin", type=float, default=None, + help="Applying SNR threshold other than 4 (full sample) or 5.5 (cosmo sample). Make sure that corresponding halo catalogue is ready" +) +parser.add_argument( + "--high-accuracy", action="store_true", help="If using a high accuracy lensed map, choose this option. (available for websky and sehgal)" ) parser.add_argument( - "--high-accuracy", action="store_true", help="If using a high accuracy websky map, choose this option." + "--inpaint", type=float, default=None, help="Inpainting for gradient leg. Enter the inpainting radius in arcmin." ) args = parser.parse_args() if args.which_sim == "websky": output_path = paths.websky_output_path elif args.which_sim == "sehgal": output_path = paths.sehgal_output_path elif args.which_sim == "agora": output_path = paths.agora_output_path -# output_path = paths.{args.which_sim}_output_path save_name = args.save_name save_dir = f"{output_path}/{save_name}" @@ -63,55 +69,112 @@ simsuite_path = f"{paths.simsuite_path}/{args.which_sim}_{paths.simsuite_version}/" cat_path = f"{paths.cat_path}/{paths.nemosim_version}/" print(" ::: catalogue:", args.which_sim, args.which_cat, "[ simsuite ver:", paths.simsuite_version, "/ nemosim ver:", paths.nemosim_version, "]") -print(" ::: hres:", args.hres_choice, "/ grad:", args.grad_choice) + +if args.cmb_tsz: print(" ::: this is CMB + tSZ run") +elif args.cmb_msub: print(" ::: this is BASELINE run") +else: print(" ::: this is CMB ONLY run") +# for now without cmb_tsz or cmb_msub, it will be cmb_only for hres as a default if args.is_meanfield: print(" ::: this is a mean-field run") # SIM SETTING ------------------------------------------------------------------ -xlmin = 200 -xlmax = 2000 - -ylmin = 200 -ylmax = 3500 # low lmax cut (high lmax cut is 6000) -ylcut = 2 - -klmin = 200 -klmax = 3000 # low lmax cut (high lmax cut is 5000) -print(" ::: hlmax =", ylmax, "and klmax =", klmax) - -tap_per = 12.0 +# for grad leg +xlmin = 200 +xlmax = 2000 + +# for hres leg +ylmin = 200 +ylmax = 3500 # low lmax cut (high lmax cut is 6000) +lycut = 2 + +# for kappa +klmin = 200 +klmax = 3000 # low lmax cut (high lmax cut is 5000) +lstep = 200 + +# for ILC +ilcmin = 200 +ilcmax = 8000 +ilc_beam_fwhm150 = 1.5 +ilc_beam_fwhm090 = 2.2 + +# for PS fitting +highres_fiducial_rms = 15 +highres_fit_ellmin = 500 +highres_fit_ellmax = 8000 +gradient_fiducial_rms = 35 +gradient_fit_ellmin = 200 +gradient_fit_ellmax = 3000 +highres_fiducial_lknee = 3000 +highres_fiducial_alpha = -4 + +# for taper +tap_per = 20.0 pad_per = 3.0 + +# for beam fwhm = 1.5 -nlevel = 15.0 +fwhm_plc = 5.0 +# for white noise +nlevel = 15.0 +nlevel_plc = 35.0 + +# for stamp cutout px = 0.5 width_deg = 120./60. maxr = width_deg * utils.degree / 2.0 # for radially binned kappa profile arcmax = 15.0 -arcstep = 1.5 +arcstep = 1.5 bin_edges = np.arange(0, arcmax, arcstep) centers = (bin_edges[1:] + bin_edges[:-1]) / 2.0 +# for kappa L +ell_edges = np.arange(klmin, klmax, lstep) +ell_cents = (ell_edges[1:] + ell_edges[:-1]) / 2.0 + +print(" ::: grad lmin =", xlmin, "and grad lmax =", xlmax) +print(" ::: hres lmin =", ylmin, "and hres lmax =", ylmax) +print(" ::: kappa Lmin =", klmin, "and kappa Lmax =", klmax) +print(" ::: theta bin size =", arcstep, "arcmin") +print(" ::: L bin size =", lstep) + def bin(data, modrmap, bin_edges): binner = stats.bin2D(modrmap, bin_edges) cents, ret = binner.bin(data) return ret +if args.inpaint is not None: + inpaint_radius = args.inpaint # default is 4 arcmin + print(" ::: inpainting for gradient leg with a radius of", inpaint_radius, "arcmin") # READING CATALOGUE ------------------------------------------------------------ if args.which_cat == "halo": - if not args.full_sample: - cat = cat_path + f"{args.which_sim}_halo_snr5p5.txt" + if args.full_sample: + cat = cat_path + f"{args.which_sim}_halo_true.txt" + elif args.highsnr_sample: + cat = cat_path + f"{args.which_sim}_halo_snr7_true.txt" else: - cat = cat_path + f"{args.which_sim}_halo.txt" + cat = cat_path + f"{args.which_sim}_halo_snr5p5_true.txt" # matching true tsz cat below + ras, decs, zs, masses = np.loadtxt(cat, unpack=True) +elif args.which_cat == "tsz_true": + if args.full_sample: + cat = cat_path + f"{args.which_sim}_tsz_true.txt" + elif args.highsnr_sample: + cat = cat_path + f"{args.which_sim}_tsz_snr7_true.txt" + else: + cat = cat_path + f"{args.which_sim}_tsz_snr5p5_true.txt" # M200c, true mass, lensed coordinates for Agora + + ras, decs, zs, masses = np.loadtxt(cat, unpack=True) + elif args.which_cat == "tsz": if args.which_sim == "websky": cat = paths.websky_tsz_cat @@ -119,42 +182,52 @@ def bin(data, modrmap, bin_edges): ras = hdu[1].data["RADeg"] decs = hdu[1].data["DECDeg"] masses = hdu[1].data["M200m"] # 1e14 Msun: websky halo cat only provides M200m - snr = hdu[1].data["SNR"] + snr = hdu[1].data["SNR"] # fixed_SNR zs = hdu[1].data["redshift"] - elif args.which_sim == "sehgal" or args.which_sim == "agora": - cat = paths.sehgal_tsz_cat + else: + if args.which_sim == "sehgal": cat = paths.sehgal_tsz_cat + elif args.which_sim == "agora": cat = paths.agora_tsz_cat + hdu = fits.open(cat) ras = hdu[1].data["RADeg"] decs = hdu[1].data["DECDeg"] - masses = hdu[1].data["M200c"] # 1e14 Msun - snr = hdu[1].data["SNR"] + masses = hdu[1].data["M200c"] # 1e14 Msun #FIXME: make it as an option + # masses = hdu[1].data["M200m"] # 1e14 Msun ##### + snr = hdu[1].data["SNR"] # fixed_SNR zs = hdu[1].data["redshift"] # SNR cut corresponding to cosmo sample is 5.5 if not args.full_sample: - snr_cut = 5.5 - keep = snr > snr_cut + if args.snmin is not None: + snr_cut = args.snmin + else: + snr_cut = 5.5 + + keep = (snr > snr_cut) + ras = ras[keep] decs = decs[keep] masses = masses[keep] zs = zs[keep] snr = snr[keep] - print(" ::: min and max SNR = %.2f and %.2f" %(snr.min(), snr.max())) -print(" ::: min and max redshift = %.2f and %.2f" %(zs.min(), zs.max()), "(mean = %.2f)" %zs.mean()) -print(" ::: min and max M200 = %.2f and %.2f" %(masses.min(), masses.max()), "(mean = %.2f)" %masses.mean()) + print(f" ::: min and max SNR = {snr.min():.2f} and {snr.max():.2f} (SNR = {snr.mean():.2f} \xB1 {snr.std():.2f})") +print(f" ::: min and max redshift = {zs.min():.2f} and {zs.max():.2f} (z = {zs.mean():.2f} \xB1 {zs.std():.2f})") +print(f" ::: min and max M200c = {masses.min():.2f} and {masses.max():.2f} (M = {masses.mean():.2f} \xB1 {masses.std():.2f})") + if args.is_meanfield: Nx = 100 * len(ras) + if args.which_cat == "tsz_true": args.which_cat = "tsz" # load random catalogue - created by mapcat.py + randcat.py cat = cat_path + f"{args.which_sim}_{args.which_cat}_randoms.txt" ras, decs = np.loadtxt(cat, unpack=True) - if not args.full_sample: - ras = ras[:Nx] - decs = decs[:Nx] + # if not args.full_sample: + ras = ras[:Nx] + decs = decs[:Nx] print(" ::: reading random catalogue for mean-field run") @@ -162,8 +235,6 @@ def bin(data, modrmap, bin_edges): print(" ::: total number of clusters for stacking =", len(ras)) - - # READING MAPS ----------------------------------------------------------------- if args.which_sim == "websky": @@ -172,52 +243,152 @@ def bin(data, modrmap, bin_edges): else: true = paths.websky_sim_path + paths.websky_kappa_reproj elif args.which_sim == "sehgal": - true = paths.sehgal_sim_path + paths.sehgal_kappa_reproj + if args.high_accuracy: + true = paths.sehgal_sim_path + paths.sehgal_dkappa_reproj + else: + true = paths.sehgal_sim_path + paths.sehgal_kappa_reproj elif args.which_sim == "agora": true = paths.agora_sim_path + paths.agora_kappa_reproj -if not args.is_observed: - print(" ::: preparing for SIGNAL maps") - - cmb = f"{simsuite_path}scmb.fits" - cmb_cib = f"{simsuite_path}scmb_cib.fits" - cmb_ksz_cib = f"{simsuite_path}scmb_ksz_cib.fits" - - cmb_tsz = f"{simsuite_path}scmb_tsz.fits" - cmb_tsz_cib = f"{simsuite_path}scmb_tsz_cib.fits" - cmb_tsz_ksz_cib = f"{simsuite_path}scmb_tsz_ksz_cib.fits" -else: - print(" ::: preparing for OBSERVED maps") - cmb = f"{simsuite_path}ocmb.fits" - cmb_cib = f"{simsuite_path}ocmb_cib.fits" - cmb_ksz_cib = f"{simsuite_path}ocmb_ksz_cib.fits" +print(" ::: preparing for OBSERVED maps") - cmb_tsz = f"{simsuite_path}ocmb_tsz.fits" - cmb_tsz_cib = f"{simsuite_path}ocmb_tsz_cib.fits" - cmb_tsz_ksz_cib = f"{simsuite_path}ocmb_tsz_ksz_cib.fits" +h_cmb = f"{simsuite_path}h_ocmb.fits" +g_cmb = f"{simsuite_path}g_ocmb.fits" print(" ::: reading true kappa map:", true) -print(" ::: reading lensed cmb map:", cmb) -print(" ::: reading lensed cmb + cib map:", cmb_cib) -print(" ::: reading lensed cmb + ksz + cib map:", cmb_ksz_cib) -print(" ::: reading lensed cmb + tsz map:", cmb_tsz) -print(" ::: reading lensed cmb + tsz + cib map:", cmb_tsz_cib) -print(" ::: reading lensed cmb + tsz + ksz + cib map:", cmb_tsz_ksz_cib) +print(" ::: reading lensed cmb map for hres:", h_cmb) +print(" ::: reading lensed cmb map for grad:", g_cmb) true_map = enmap.read_map(true, delayed=False) -cmb_map = enmap.read_map(cmb, delayed=False) -cmb_cib_map = enmap.read_map(cmb_cib, delayed=False) -cmb_ksz_cib_map = enmap.read_map(cmb_ksz_cib, delayed=False) -cmb_tsz_map = enmap.read_map(cmb_tsz, delayed=False) -cmb_tsz_cib_map = enmap.read_map(cmb_tsz_cib, delayed=False) -cmb_tsz_ksz_cib_map = enmap.read_map(cmb_tsz_ksz_cib, delayed=False) +h_cmb_map = enmap.read_map(h_cmb, delayed=False) +g_cmb_map = enmap.read_map(g_cmb, delayed=False) + +if args.cmb_tsz: + + cmb_tsz150 = f"{simsuite_path}ocmb_tsz150.fits" + cmb_tsz090 = f"{simsuite_path}ocmb_tsz090.fits" + + print(" ::: reading lensed cmb + tsz map at 150GHz:", cmb_tsz150) + print(" ::: reading lensed cmb + tsz map at 90GHz:", cmb_tsz090) + + cmb_tsz_map150 = enmap.read_map(cmb_tsz150, delayed=False) + cmb_tsz_map090 = enmap.read_map(cmb_tsz090, delayed=False) + +if args.cmb_msub: + + cmb_tsz090 = f"{simsuite_path}ocmb_tsz090_actfoot.fits" + cmb_tsz150 = f"{simsuite_path}ocmb_tsz150_actfoot.fits" + + cmb_tsz_map090 = enmap.read_map(cmb_tsz090, delayed=False) + cmb_tsz_map150 = enmap.read_map(cmb_tsz150, delayed=False) + + model_path = f"/home3/eunseong/nemo-sims/nemo-sim-kit48/agora_cmb_tsz/" + + model090 = model_path + f"clusterModelMaps/clusterModelMap_f090.fits" + model150 = model_path + f"clusterModelMaps/clusterModelMap_f150.fits" + + print(" ::: reading model image map:", model090) + print(" ::: reading model image map:", model150) + + model_map090 = enmap.read_map(model090, delayed=False) + model_map150 = enmap.read_map(model150, delayed=False) + + print(" ::: model image subtraction!") + + cmb_msub_map090 = cmb_tsz_map090 - model_map090 + cmb_msub_map150 = cmb_tsz_map150 - model_map150 print(" ::: maps are ready!") +# FOR ILC / COADD ---------------------------------------------------------------- + +def fit_p1d( + l_edges, cents, p1d, which, xout, bfunc1, bfunc2, rms=None, lmin=None, lmax=None +): + # function for fitting 1D power spectrum of given stamp + b1 = bfunc1 if bfunc1 is not None else lambda x: 1 + b2 = bfunc2 if bfunc2 is not None else lambda x: 1 + + tfunc = lambda x: theory.lCl("TT", x) * b1(x) * b2(x) + + # PS fitting + # Select region for fit + sel = np.logical_and(cents > lmin, cents < lmax) + delta_ells = np.diff(l_edges)[sel] + ells = cents[sel] + cls = p1d[sel] + cltt = tfunc(ells) # fiducial Cltt + + if which == "plc" or which == "act" or which == "act_cross": + w0 = gradient_fiducial_rms if which=='plc' else highres_fiducial_rms + sigma2 = stats.get_sigma2(ells, cltt, w0, delta_ells, fsky, ell0=0, alpha=1) + func = stats.fit_cltt_power(ells, cls, tfunc, w0, sigma2, ell0=0, alpha=1, fix_knee=True) + + elif which == "apcross": + w0 = gradient_fiducial_rms + w0p = highres_fiducial_rms + sigma2 = stats.get_sigma2(ells, cltt, w0, delta_ells, fsky, ell0=0, alpha=0, w0p=w0p, ell0p=0, alphap=1, clxx=cltt, clyy=cltt) + func = stats.fit_cltt_power(ells, cls, tfunc, w0, sigma2, ell0=0, alpha=1, fix_knee=True) + + ret = func(xout) + + ret[xout < 2] = 0 + assert np.all(np.isfinite(ret)) + return ret + + +def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): + # A simple two array ILC solution + # Returns beam deconvolved + sel = np.logical_and(modlmap >= ilcmin, modlmap <= ilcmax) + nells = modlmap[sel].size + cov = np.zeros((nells, 2, 2)) + cov[:, 0, 0] = p11[sel] + cov[:, 1, 1] = p22[sel] + cov[:, 0, 1] = p12[sel] + cov[:, 1, 0] = p12[sel] + ms = np.stack([m1[sel], m2[sel]]).swapaxes(0, 1) + rs = np.stack([b1[sel], b2[sel]]).swapaxes(0, 1) + num = np.linalg.solve(cov, ms) + den = np.linalg.solve(cov, rs) + tcov = 1.0 / np.einsum("ij,ij->i", rs, den) + ksolve = np.einsum("ij,ij->i", rs, num) * tcov + assert np.all(np.isfinite(ksolve)) + ret = m1 * 0 + ret[sel] = ksolve + tret = p11 * 0 + tret[sel] = tcov + return ret, tret + + +##### TEST ##### --------------------------------------------------------------- + +rstamp = 120.8 * utils.arcmin # to match the stamp size from the sim (241, 241) +res = 0.5 * utils.arcmin + +lmin = 200 +lmax = 3000 # 20000 +rmin = 0. +rmax = 10 * utils.arcmin +rwidth = 1.0 * utils.arcmin + +# agora +M200c = 6.37e14 +z = 0.6 + +# this returns conc = 3.326 using Klypin et al 2016 c(M,z) +# k2dmap is our fixed template on 2D map - this is tapered and filtered +_,_,_,_,_,_,_,_,k2dmap,_ = lensing.kappa_nfw_profiley(mass=M200c,conc=None,z=z,z_s=1100.,background='critical',delta=200,apply_filter=True,lmin=lmin,lmax=lmax,res=res,rstamp=rstamp,rmin=rmin,rmax=rmax,rwidth=rwidth) + +##### TEST ##### --------------------------------------------------------------- + + + + # LOOP OVER ASSIGNED TASKS ----------------------------------------------------- if args.is_test: @@ -233,7 +404,8 @@ def bin(data, modrmap, bin_edges): i = task cper = int((j + 1) / len(my_tasks) * 100.0) if rank==0: print(f"Rank {rank} performing task {task} as index {j} ({cper}% complete.).") - + + coords = np.array([decs[i], ras[i]]) * utils.degree # cut out a stamp from the simulated map @@ -247,65 +419,68 @@ def bin(data, modrmap, bin_edges): pixwin=False ) # true kappa - cmb = reproject.thumbnails( - cmb_map, + h_cmb = reproject.thumbnails( + h_cmb_map, coords, r=maxr, res=px * utils.arcmin, proj="tan", oversample=2, pixwin=True if args.high_accuracy else False # only True for maps made natively in CAR - ) # lensed cmb + ) # lensed cmb for hres - cmb_cib = reproject.thumbnails( - cmb_cib_map, + g_cmb = reproject.thumbnails( + g_cmb_map, coords, r=maxr, res=px * utils.arcmin, proj="tan", oversample=2, - pixwin=False - ) # lensed cmb + cib - - cmb_ksz_cib = reproject.thumbnails( - cmb_ksz_cib_map, - coords, - r=maxr, - res=px * utils.arcmin, - proj="tan", - oversample=2, - pixwin=False - ) # lensed cmb + ksz + cib - - cmb_tsz = reproject.thumbnails( - cmb_tsz_map, - coords, - r=maxr, - res=px * utils.arcmin, - proj="tan", - oversample=2, - pixwin=False - ) # lensed cmb + tsz - - cmb_tsz_cib = reproject.thumbnails( - cmb_tsz_cib_map, - coords, - r=maxr, - res=px * utils.arcmin, - proj="tan", - oversample=2, - pixwin=False - ) # lensed cmb + tsz + cib + pixwin=True if args.high_accuracy else False # only True for maps made natively in CAR + ) # lensed cmb for grad + + if args.cmb_tsz: + cmb_tsz150 = reproject.thumbnails( + cmb_tsz_map150, + coords, + r=maxr, + res=px * utils.arcmin, + proj="tan", + oversample=2, + pixwin=False + ) # lensed cmb + tsz150 + + cmb_tsz090 = reproject.thumbnails( + cmb_tsz_map090, + coords, + r=maxr, + res=px * utils.arcmin, + proj="tan", + oversample=2, + pixwin=False + ) # lensed cmb + tsz090 + + if args.cmb_msub: + cmb_msub150 = reproject.thumbnails( + cmb_msub_map150, + coords, + r=maxr, + res=px * utils.arcmin, + proj="tan", + oversample=2, + pixwin=False + ) # lensed cmb + tsz150 - model150 + + cmb_msub090 = reproject.thumbnails( + cmb_msub_map090, + coords, + r=maxr, + res=px * utils.arcmin, + proj="tan", + oversample=2, + pixwin=False + ) # lensed cmb + tsz090 - model090 - cmb_tsz_ksz_cib = reproject.thumbnails( - cmb_tsz_ksz_cib_map, - coords, - r=maxr, - res=px * utils.arcmin, - proj="tan", - oversample=2, - pixwin=False - ) # lensed cmb + tsz + ksz + cib # initialise calculations based on geometry if j == 0: @@ -315,9 +490,17 @@ def bin(data, modrmap, bin_edges): wcs = kstamp.wcs modlmap = enmap.modlmap(shape, wcs) modrmap = enmap.modrmap(shape, wcs) - - assert wcsutils.equal(kstamp.wcs, cmb.wcs) - assert wcsutils.equal(kstamp.wcs, cmb_tsz.wcs) + + assert wcsutils.equal(kstamp.wcs, h_cmb.wcs) + assert wcsutils.equal(kstamp.wcs, g_cmb.wcs) + + if args.cmb_tsz: + assert wcsutils.equal(kstamp.wcs, cmb_tsz150.wcs) + assert wcsutils.equal(kstamp.wcs, cmb_tsz090.wcs) + + if args.cmb_msub: + assert wcsutils.equal(kstamp.wcs, cmb_msub150.wcs) + assert wcsutils.equal(kstamp.wcs, cmb_msub090.wcs) # get an edge taper map and apodize taper = maps.get_taper( @@ -330,61 +513,321 @@ def bin(data, modrmap, bin_edges): taper = taper[0] # evaluate the 2D Gaussian beam on an isotropic Fourier grid - beam2d = maps.gauss_beam(modlmap, fwhm) + beam2d = maps.gauss_beam(modlmap, fwhm) + beam2d_plc = maps.gauss_beam(modlmap, fwhm_plc) + + if args.cmb_tsz or args.cmb_msub: + bfunc150 = lambda x: maps.gauss_beam(ilc_beam_fwhm150, x) + bfunc90 = lambda x: maps.gauss_beam(ilc_beam_fwhm090, x) + + act_150_kbeam2d = bfunc150(modlmap) + act_90_kbeam2d = bfunc90(modlmap) + plc_kbeam2d = beam2d_plc # build Fourier space masks for lensing reconstruction xmask = maps.mask_kspace(shape, wcs, lmin=xlmin, lmax=xlmax) # grad - ymask = maps.mask_kspace(shape, wcs, lmin=ylmin, lmax=ylmax, lxcut=ylcut, lycut=ylcut) # hres + ymask = maps.mask_kspace(shape, wcs, lmin=ylmin, lmax=ylmax, lxcut=lycut, lycut=lycut) # hres kmask = maps.mask_kspace(shape, wcs, lmin=klmin, lmax=klmax) # kappa # get theory spectrum and build interpolated 2D Fourier CMB from theory and maps theory = cosmology.default_theory() ucltt2d = theory.lCl("TT", modlmap) - # total spectrum includes beam-deconvolved noise - npower = (nlevel * np.pi/180./60.)**2. - tcltt2d = ucltt2d + npower/beam2d**2. - tcltt2d[~np.isfinite(tcltt2d)] = 0 + if args.cmb_tsz or args.cmb_msub: + # bin size and range for 1D binned power spectrum + minell = 2 * maps.minimum_ell(shape, wcs) + l_edges = np.arange(minell / 2, 8001, minell) + lbinner = stats.bin2D(modlmap, l_edges) + # PS correction factor + w2 = np.mean(taper ** 2) + # fsky for bandpower variance + fsky = enmap.area(shape, wcs) * w2 / 4.0 / np.pi + - # same filter as the post-reconstuction for true kappa - k_stamp = maps.filter_map(kstamp, kmask) + # same filter as the post-reconstuction for true kappa (also tapered!) + k_stamp = maps.filter_map(kstamp*taper, kmask) s.add_to_stack("kstamp", k_stamp) binned_true = bin(k_stamp, modrmap * (180 * 60 / np.pi), bin_edges) s.add_to_stats("tk1d", binned_true) - # choose the map for each leg - hres = globals()[args.hres_choice] - grad = globals()[args.grad_choice] - # taper stamp - tapered_hres = hres * taper - tapered_grad = grad * taper + if args.cmb_tsz: + hres150 = cmb_tsz150 + hres90 = cmb_tsz090 + elif args.cmb_msub: + hres150 = cmb_msub150 + hres90 = cmb_msub090 + else: + hres = h_cmb + + + # to filter weird stamps when model image subtraction is done for the mean-field + if args.cmb_msub: + if np.any(abs(hres150) > 1e3) or np.any(abs(hres90) > 1e3): + print(f"{task} has anomalously high hres150 or hres90") + continue - if not args.is_observed: - # get a Fourier transformed stamp - k_hres = enmap.fft(tapered_hres, normalize="phys") + # taper stamp + tapered_grad = g_cmb * taper + + if args.cmb_tsz or args.cmb_msub: + tapered_hres150 = hres150 * taper + tapered_hres90 = hres90 * taper + else: + tapered_hres = hres * taper + + s.add_to_stack("grad2d_before", g_cmb) + g_filtered = maps.filter_map(g_cmb, ymask) + s.add_to_stack("grad2d_before_filtered", g_filtered) + binned_grad = bin(g_cmb, g_cmb.modrmap() * (180 * 60 / np.pi), bin_edges) + s.add_to_stats("grad1d_before", binned_grad) + + + # this is the old inpainting methood + if args.inpaint: # for gradient leg (should be used with args.is_observed at the moment) + """ + If inpainting, we + (1) resample the stamp to 64x64 (2 arcmin pixels) + (2) Inpaint a hole of radius 4 arcmin + """ + rmin = inpaint_radius * utils.arcmin + crop_pixels = int(inpaint_radius * 4. / px) # 16 arcminutes wide for 4 arcmin radius + cutout = maps.crop_center(tapered_grad, cropy=crop_pixels, cropx=crop_pixels, sel=False) + cutout_sel = maps.crop_center(tapered_grad, cropy=crop_pixels, cropx=crop_pixels, sel=True) + Ndown, Ndown2 = cutout.shape[-2:] + if Ndown != Ndown2: raise Exception + + if j==0: # set up for planck-like beam and noise + from orphics import pixcov + pshape = cutout.shape + pwcs = cutout.wcs + beam_fn = lambda x: maps.gauss_beam(fwhm_plc, x) + ipsizemap = enmap.pixsizemap(pshape, pwcs) + pivar = maps.ivar(pshape, pwcs, nlevel_plc, ipsizemap=ipsizemap) + pcov = pixcov.tpcov_from_ivar(Ndown, pivar, theory.lCl, beam_fn) + geo = pixcov.make_geometry(pshape, pwcs, rmin, n=Ndown, deproject=True, iau=False, res=None, pcov=pcov) + + cutout = pixcov.inpaint_stamp(cutout, geo) + tapered_grad[cutout_sel] = cutout.copy() + tapered_grad = np.nan_to_num(tapered_grad) + g_cmb = tapered_grad / taper + g_cmb = np.nan_to_num(g_cmb) + + + s.add_to_stack("grad2d_after", g_cmb) + g_filtered = maps.filter_map(g_cmb, ymask) + s.add_to_stack("grad2d_after_filtered", g_filtered) + binned_grad = bin(g_cmb, g_cmb.modrmap() * (180 * 60 / np.pi), bin_edges) + s.add_to_stats("grad1d_after", binned_grad) + + + + if args.cmb_tsz or args.cmb_msub: + + # ILC will return a beam-deconvolved Fourier transformed stamp k_grad = enmap.fft(tapered_grad, normalize="phys") + k150 = enmap.fft(tapered_hres150, normalize="phys") + k90 = enmap.fft(tapered_hres90, normalize="phys") + assert np.all(np.isfinite(k150)) + assert np.all(np.isfinite(k90)) + + # Fourier map -> PS + pow = lambda x, y: (x * y.conj()).real + + # measure the binned power spectrum from given stamp + act_cents, act_p1d_150 = lbinner.bin(pow(k150, k150) / w2) + act_cents, act_p1d_90 = lbinner.bin(pow(k90, k90) / w2) + act_cents, act_p1d_150_90 = lbinner.bin(pow(k150, k90) / w2) + plc_cents, plc_p1d = lbinner.bin(pow(k_grad, k_grad) / w2) + + # fit power spectra + tclaa_150 = fit_p1d( + l_edges, + act_cents, + act_p1d_150, + "act", + modlmap, + bfunc150, + bfunc150, + rms=highres_fiducial_rms, + lmin=highres_fit_ellmin, + lmax=highres_fit_ellmax, + ) + + tclaa_90 = fit_p1d( + l_edges, + act_cents, + act_p1d_90, + "act", + modlmap, + bfunc90, + bfunc90, + rms=highres_fiducial_rms, + lmin=highres_fit_ellmin, + lmax=highres_fit_ellmax, + ) + + tclaa_150_90 = fit_p1d( + l_edges, + act_cents, + act_p1d_150_90, + "act_cross", + modlmap, + bfunc150, + bfunc90, + rms=0, + lmin=highres_fit_ellmin, + lmax=highres_fit_ellmax, + ) + + tclpp = fit_p1d( + l_edges, + plc_cents, + plc_p1d, + "plc", + modlmap, + lambda x: maps.gauss_beam(x, fwhm_plc), + lambda x: maps.gauss_beam(x, fwhm_plc), + rms=gradient_fiducial_rms, + lmin=gradient_fit_ellmin, + lmax=gradient_fit_ellmax, + ) + + # ILC / coadd + k_hres, tclaa = ilc( + modlmap, + k150, + k90, + tclaa_150, + tclaa_90, + tclaa_150_90, + act_150_kbeam2d, + act_90_kbeam2d, + ) # beam deconvolved + + k_grad = k_grad / plc_kbeam2d + tclpp = tclpp / (plc_kbeam2d ** 2.0) + + # # fit cross-power of gradient and ILC hres + # cents, c_ap = lbinner.bin(pow(k_hres, k_grad) / w2) + # tclap = fit_p1d( + # l_edges, + # cents, + # c_ap, + # "apcross", + # modlmap, + # None, + # None, + # rms=0, + # lmin=highres_fit_ellmin, + # lmax=gradient_fit_ellmax + # ) # shouldn't this be beam deconvolved as well? + # tclap[~np.isfinite(tclap)] = 0 + else: + # cmb only hres run doesn't involve ILC + # total spectrum includes beam-deconvolved noise + npower = (nlevel * np.pi/180./60.)**2. + npower_plc = (nlevel_plc * np.pi/180./60.)**2. + + tclaa = ucltt2d + npower/beam2d**2. + tclpp = ucltt2d + npower_plc/beam2d_plc**2. + # get a beam-deconvolved Fourier transformed stamp + k_grad = enmap.fft(tapered_grad, normalize="phys")/beam2d_plc k_hres = enmap.fft(tapered_hres, normalize="phys")/beam2d - k_grad = enmap.fft(tapered_grad, normalize="phys")/beam2d - assert np.all(np.isfinite(k_hres)) assert np.all(np.isfinite(k_grad)) + assert np.all(np.isfinite(k_hres)) + + tclaa[~np.isfinite(tclaa)] = 0 + tclpp[~np.isfinite(tclpp)] = 0 # build symlens dictionary feed_dict = { "uC_T_T" : ucltt2d, - "tC_T_T" : tcltt2d, + "tC_A_T_A_T": tclaa, # approximate ACT power spectrum, ACT beam deconvolved + "tC_P_T_P_T": tclpp, # approximate Planck power spectrum, Planck beam deconvolved + "tC_A_T_P_T": ucltt2d, # same lensed theory as above, no instrumental noise + "tC_P_T_A_T": ucltt2d, # same lensed theory as above, no instrumental noise "X" : k_grad, # grad leg "Y" : k_hres, # hres leg } - # do lensing reconstruction in Fourier space - rkmap = qe.reconstruct(shape, wcs, feed_dict, estimator="hdv", XY="TT", xmask=xmask, ymask=ymask, kmask=kmask, physical_units=True) - + # Sanity check + for key in feed_dict.keys(): + assert np.all(np.isfinite(feed_dict[key])) + + # ask for reconstruction in Fourier space + cqe = symlens.QE( + shape, + wcs, + feed_dict, + estimator="hdv", + XY="TT", + xmask=xmask, + ymask=ymask, + field_names=["P", "A"], + groups=None, + kmask=kmask, + ) + + # Fourier space lens reconstruction + rkmap = cqe.reconstruct(feed_dict, xname="X_l1", yname="Y_l2", physical_units=True) + assert np.all(np.isfinite(rkmap)) - + + + # for cross-correlation check ---------------------------------------------- + + ikmap = enmap.fft(k_stamp, normalize="phys") # k_stamp is tapered and filtered true kappa + assert np.all(np.isfinite(ikmap)) + + p2d_auto = (ikmap * ikmap.conj()).real + p2d_cros = (ikmap * rkmap.conj()).real + + p1d_auto = bin(p2d_auto, modlmap, ell_edges) + p1d_cros = bin(p2d_cros, modlmap, ell_edges) + + s.add_to_stats("p1d_auto", p1d_auto) + s.add_to_stats("p1d_cros", p1d_cros) + s.add_to_stack("true_real", ikmap.real) + s.add_to_stack("true_imag", ikmap.imag) + s.add_to_stack("recon_real", rkmap.real) + s.add_to_stack("recon_imag", rkmap.imag) + + # for cross-correlation check ---------------------------------------------- + + ktemp = k2dmap # this is tapered and filtered template (in real space) + + template = enmap.fft(ktemp, normalize="phys") + assert np.all(np.isfinite(template)) + + p2d_cross_t = (ikmap * template.conj()).real # cross with true kappa + p2d_cross_r = (rkmap * template.conj()).real # cross with reconstruction - what about the meanfield here? post-processing! + p2d_auto_template = (template * template.conj()).real + + p1d_cross_t = bin(p2d_cross_t, modlmap, ell_edges) + p1d_cross_r = bin(p2d_cross_r, modlmap, ell_edges) + p1d_auto_template = bin(p2d_auto_template, modlmap, ell_edges) + + s.add_to_stats("p1d_cross_t", p1d_cross_t) + s.add_to_stats("p1d_cross_r", p1d_cross_r) + s.add_to_stats("p1d_auto_template", p1d_auto_template) + + l = ell_cents * ell_cents + s.add_to_stats("p1d_cross_t_l2", l * p1d_cross_t) + s.add_to_stats("p1d_cross_r_l2", l * p1d_cross_r) + s.add_to_stats("p1d_auto_template_l2", l * p1d_auto_template) + + s.add_to_stack("temp", k2dmap) + binned_temp = bin(k2dmap, modrmap * (180 * 60 / np.pi), bin_edges) + s.add_to_stats("temp1d", binned_temp) + + + + + # transform to real space kappa = enmap.ifft(rkmap, normalize="phys").real @@ -394,10 +837,20 @@ def bin(data, modrmap, bin_edges): s.add_to_stats("k1d", binned_kappa) # stack stamps pre-reconstruction as well (same filter as hres leg) + if args.cmb_tsz: + hres = cmb_tsz150 + elif args.cmb_msub: + hres = cmb_msub150 + else: + hres = h_cmb + grad = g_cmb + hres_stamp = maps.filter_map(hres, ymask) grad_stamp = maps.filter_map(grad, ymask) s.add_to_stack("hres_stamp_st", hres_stamp) s.add_to_stack("grad_stamp_st", grad_stamp) + s.add_to_stack("hres_noF_stamp_st", hres) + s.add_to_stack("grad_noF_stamp_st", grad) # save the list of masses and redshifts for matched stack mass fitting if not args.is_meanfield: @@ -420,16 +873,22 @@ def bin(data, modrmap, bin_edges): # stacks before lensing reconstruction hres_st = s.stacks["hres_stamp_st"] - grad_st = s.stacks["grad_stamp_st"] + grad_st = s.stacks["grad_stamp_st"] + hres_noF_st = s.stacks["hres_noF_stamp_st"] + grad_noF_st = s.stacks["grad_noF_stamp_st"] hres_zoom = hres_st[100:140,100:140] grad_zoom = grad_st[100:140,100:140] + hres_noF_zoom = hres_noF_st[100:140,100:140] + grad_noF_zoom = grad_noF_st[100:140,100:140] modrmap = hres_zoom.modrmap() modrmap = np.rad2deg(modrmap)*60. hbinned = bin(hres_zoom, modrmap, bin_edges) gbinned = bin(grad_zoom, modrmap, bin_edges) + hbinned_noF = bin(hres_noF_zoom, modrmap, bin_edges) + gbinned_noF = bin(grad_noF_zoom, modrmap, bin_edges) io.plot_img(hres_st, f"{save_dir}/{save_name}_0hres.png") io.plot_img(grad_st, f"{save_dir}/{save_name}_0grad.png") @@ -437,9 +896,19 @@ def bin(data, modrmap, bin_edges): io.plot_img(grad_zoom, f"{save_dir}/{save_name}_0grad_zoom.png") np.save(f"{save_dir}/{save_name}_0hres.npy", hres_st) np.save(f"{save_dir}/{save_name}_0grad.npy", grad_st) + io.save_cols(f"{save_dir}/{save_name}_0binned_hres.txt", (centers, hbinned)) io.save_cols(f"{save_dir}/{save_name}_0binned_grad.txt", (centers, gbinned)) + io.plot_img(hres_noF_st, f"{save_dir}/{save_name}_0hres_noF.png") + io.plot_img(grad_noF_st, f"{save_dir}/{save_name}_0grad_noF.png") + io.plot_img(hres_noF_zoom, f"{save_dir}/{save_name}_0hres_noF_zoom.png") + io.plot_img(grad_noF_zoom, f"{save_dir}/{save_name}_0grad_noF_zoom.png") + np.save(f"{save_dir}/{save_name}_0hres_noF.npy", hres_noF_st) + np.save(f"{save_dir}/{save_name}_0grad_noF.npy", grad_noF_st) + io.save_cols(f"{save_dir}/{save_name}_0binned_hres_noF.txt", (centers, hbinned_noF)) + io.save_cols(f"{save_dir}/{save_name}_0binned_grad_noF.txt", (centers, gbinned_noF)) + # reconstructed lensing field kmap = s.stacks["kstamp"] lmap = s.stacks["lstamp"] @@ -486,6 +955,20 @@ def bin(data, modrmap, bin_edges): if not args.is_meanfield: np.savetxt(f"{save_dir}/{save_name}_z_mass1e14.txt", np.c_[s.vectors["redshift"], s.vectors["masses"]]) + np.savetxt(f"{save_dir}/{save_name}_tkappa_1d_ind.txt", np.c_[s.vectors["tk1d"]]) + np.savetxt(f"{save_dir}/{save_name}_rkappa_1d_ind.txt", np.c_[s.vectors["k1d"]]) + + # for cross-correlation check ---------------------------------------------- + + if args.is_meanfield: + save_dir = save_dir + "_mf" + os.makedirs(save_dir) + + s.dump(f"{save_dir}") + + # for cross-correlation check ---------------------------------------------- + + elapsed = t.time() - start print("\r ::: entire run took %.1f seconds" %elapsed) From 570088f4c0e596f3888f8b17c7f18aa674078cc7 Mon Sep 17 00:00:00 2001 From: Eunseong Lee <36232273+eunseongleee@users.noreply.github.com> Date: Tue, 19 May 2026 10:24:08 -0400 Subject: [PATCH 04/14] checking --- recon_sim.py | 710 ++++++++++++++++++++++++++++++++++++++++++++++----- stack.py | 83 +++--- 2 files changed, 693 insertions(+), 100 deletions(-) diff --git a/recon_sim.py b/recon_sim.py index 8d14abf..2133f3f 100644 --- a/recon_sim.py +++ b/recon_sim.py @@ -5,7 +5,8 @@ import numpy as np from pixell import enmap, reproject, utils, wcsutils, bunch from orphics import mpi, maps, stats, io, cosmology, lensing -import symlens +# import symlens +from symlens.qe import QE import argparse import sys, os @@ -28,10 +29,40 @@ "which_cat", type=str, help="Choose the catalogue type e.g. halo or tsz." ) parser.add_argument( - "--cmb-tsz", action="store_true", help="CMB + tSZ map for high resolution leg." + "--cmb", action="store_true", help="CMB for high resolution leg." ) parser.add_argument( - "--cmb-msub", action="store_true", help="CMB + tSZ and model image subtraction for high resolution leg." + "--cmb-tsz", action="store_true", help="CMB+tSZ for high resolution leg." +) +parser.add_argument( + "--cmb-msub", action="store_true", help="CMB+tSZ-model for high resolution leg." +) +parser.add_argument( + "--cmb-ksz", action="store_true", help="CMB+kSZ for hres and CMB+kSZ for gradient leg." +) +parser.add_argument( + "--cmb-cib", action="store_true", help="CMB+CIB(90,150) for hres and CMB+CIB150 for gradient leg." +) +parser.add_argument( + "--cmb-ksz-tsz", action="store_true", help="CMB+kSZ+tSZ for hres and CMB+kSZ for gradient leg." +) +parser.add_argument( + "--cmb-cib-tsz", action="store_true", help="CMB+CIB+tSZ(90,150) for hres and CMB+CIB150 for gradient leg." +) +parser.add_argument( + "--cmb-ksz-msub", action="store_true", help="CMB+kSZ+tSZ-model for hres and CMB+kSZ for gradient leg." +) +parser.add_argument( + "--cmb-cib-msub", action="store_true", help="CMB+CIB+tSZ(90,150)-model for hres and CMB+CIB150 for gradient leg." +) +parser.add_argument( + "--cmb-ksz-cib", action="store_true", help="CMB+kSZ+CIB(90,150) for hres and CMB+kSZ+CIB150 for gradient leg." +) +parser.add_argument( + "--cmb-ksz-cib-tsz", action="store_true", help="CMB+kSZ+CIB+tSZ(90,150) for hres and CMB+kSZ+CIB150 for gradient leg." +) +parser.add_argument( + "--cmb-ksz-cib-msub", action="store_true", help="CMB+kSZ+CIB+tSZ-model(90,150) for hres and CMB+kSZ+CIB150 for gradient leg." ) parser.add_argument( "--is-meanfield", action="store_true", help="This is a mean-field run." @@ -70,10 +101,18 @@ cat_path = f"{paths.cat_path}/{paths.nemosim_version}/" print(" ::: catalogue:", args.which_sim, args.which_cat, "[ simsuite ver:", paths.simsuite_version, "/ nemosim ver:", paths.nemosim_version, "]") -if args.cmb_tsz: print(" ::: this is CMB + tSZ run") -elif args.cmb_msub: print(" ::: this is BASELINE run") -else: print(" ::: this is CMB ONLY run") -# for now without cmb_tsz or cmb_msub, it will be cmb_only for hres as a default +if args.cmb: print(" ::: this is CMB only run") +elif args.cmb_tsz: print(" ::: this is CMB + tSZ run") +elif args.cmb_msub: print(" ::: this is (CMB + tSZ - model) run") +elif args.cmb_ksz: print(" ::: this is CMB + kSZ run") +elif args.cmb_cib: print(" ::: this is CMB + CIB run") +elif args.cmb_ksz_cib: print(" ::: this is CMB + kSZ + CIB run") +elif args.cmb_ksz_tsz: print(" ::: this is CMB + kSZ + tSZ run") +elif args.cmb_cib_tsz: print(" ::: this is CMB + CIB + tSZ run") +elif args.cmb_ksz_msub: print(" ::: this is CMB + kSZ + tSZ - model run") +elif args.cmb_cib_msub: print(" ::: this is CMB + CIB + tSZ - model run") +elif args.cmb_ksz_cib_tsz: print(" ::: this is CMB + tSZ + kSZ + CIB run") +elif args.cmb_ksz_cib_msub: print(" ::: this is BASELINE (CMB + tSZ + kSZ + CIB - model) run") if args.is_meanfield: print(" ::: this is a mean-field run") @@ -111,10 +150,9 @@ highres_fiducial_alpha = -4 # for taper -tap_per = 20.0 +tap_per = 20.0 #12.0 pad_per = 3.0 - # for beam fwhm = 1.5 fwhm_plc = 5.0 @@ -130,7 +168,7 @@ # for radially binned kappa profile arcmax = 15.0 -arcstep = 1.5 +arcstep = 1.5 # default is 1.5; test with 2.25 and 3.0 bin_edges = np.arange(0, arcmax, arcstep) centers = (bin_edges[1:] + bin_edges[:-1]) / 2.0 @@ -161,6 +199,7 @@ def bin(data, modrmap, bin_edges): elif args.highsnr_sample: cat = cat_path + f"{args.which_sim}_halo_snr7_true.txt" else: + # cat = cat_path + f"{args.which_sim}_halo_snr5p5.txt" cat = cat_path + f"{args.which_sim}_halo_snr5p5_true.txt" # matching true tsz cat below ras, decs, zs, masses = np.loadtxt(cat, unpack=True) @@ -267,8 +306,8 @@ def bin(data, modrmap, bin_edges): if args.cmb_tsz: - cmb_tsz150 = f"{simsuite_path}ocmb_tsz150.fits" - cmb_tsz090 = f"{simsuite_path}ocmb_tsz090.fits" + cmb_tsz150 = f"{simsuite_path}h_ocmb_tsz150.fits" + cmb_tsz090 = f"{simsuite_path}h_ocmb_tsz090.fits" print(" ::: reading lensed cmb + tsz map at 150GHz:", cmb_tsz150) print(" ::: reading lensed cmb + tsz map at 90GHz:", cmb_tsz090) @@ -278,27 +317,147 @@ def bin(data, modrmap, bin_edges): if args.cmb_msub: - cmb_tsz090 = f"{simsuite_path}ocmb_tsz090_actfoot.fits" - cmb_tsz150 = f"{simsuite_path}ocmb_tsz150_actfoot.fits" + cmb_msub090 = f"{simsuite_path}h_ocmb_msub090_test.fits" + cmb_msub150 = f"{simsuite_path}h_ocmb_msub150_test.fits" - cmb_tsz_map090 = enmap.read_map(cmb_tsz090, delayed=False) - cmb_tsz_map150 = enmap.read_map(cmb_tsz150, delayed=False) + print(" ::: reading the model subtracted maps!", cmb_msub090) + print(" ::: reading the model subtracted maps!", cmb_msub150) + + cmb_msub_map090 = enmap.read_map(cmb_msub090, delayed=False) + cmb_msub_map150 = enmap.read_map(cmb_msub150, delayed=False) + +if args.cmb_ksz: + + h_cmb_ksz = f"{simsuite_path}h_ocmb_ksz.fits" + g_cmb_ksz = f"{simsuite_path}g_ocmb_ksz.fits" + + print(" ::: reading lensed cmb + ksz map for hres:", h_cmb_ksz) + print(" ::: reading lensed cmb + ksz map for grad:", g_cmb_ksz) + + h_cmb_ksz_map = enmap.read_map(h_cmb_ksz, delayed=False) + g_cmb_ksz_map = enmap.read_map(g_cmb_ksz, delayed=False) + +if args.cmb_cib: + + h_cmb_cib090 = f"{simsuite_path}h_ocmb_cib090_5.fits" + h_cmb_cib150 = f"{simsuite_path}h_ocmb_cib150_5.fits" + g_cmb_cib = f"{simsuite_path}g_ocmb_cib5.fits" + + print(" ::: reading lensed cmb + cib090 map for hres:", h_cmb_cib090) + print(" ::: reading lensed cmb + cib150 map for hres:", h_cmb_cib150) + print(" ::: reading lensed cmb + cib map for grad:", g_cmb_cib) + + h_cmb_cib090_map = enmap.read_map(h_cmb_cib090, delayed=False) + h_cmb_cib150_map = enmap.read_map(h_cmb_cib150, delayed=False) + g_cmb_cib_map = enmap.read_map(g_cmb_cib, delayed=False) + + +if args.cmb_ksz_tsz: + + h_cmb_ksz_tsz090 = f"{simsuite_path}h_ocmb_tsz090_ksz.fits" + h_cmb_ksz_tsz150 = f"{simsuite_path}h_ocmb_tsz150_ksz.fits" + g_cmb_ksz = f"{simsuite_path}g_ocmb_ksz.fits" - model_path = f"/home3/eunseong/nemo-sims/nemo-sim-kit48/agora_cmb_tsz/" + print(" ::: reading lensed cmb + ksz + tsz090 map for hres:", h_cmb_ksz_tsz090) + print(" ::: reading lensed cmb + ksz + tsz150 map for hres:", h_cmb_ksz_tsz150) + print(" ::: reading lensed cmb + ksz map for grad:", g_cmb_ksz) - model090 = model_path + f"clusterModelMaps/clusterModelMap_f090.fits" - model150 = model_path + f"clusterModelMaps/clusterModelMap_f150.fits" + h_cmb_ksz_tsz090_map = enmap.read_map(h_cmb_ksz_tsz090, delayed=False) + h_cmb_ksz_tsz150_map = enmap.read_map(h_cmb_ksz_tsz150, delayed=False) + g_cmb_ksz_map = enmap.read_map(g_cmb_ksz, delayed=False) - print(" ::: reading model image map:", model090) - print(" ::: reading model image map:", model150) +if args.cmb_cib_tsz: - model_map090 = enmap.read_map(model090, delayed=False) - model_map150 = enmap.read_map(model150, delayed=False) + h_cmb_cib_tsz090 = f"{simsuite_path}h_ocmb_tsz090_cib5.fits" + h_cmb_cib_tsz150 = f"{simsuite_path}h_ocmb_tsz150_cib5.fits" + g_cmb_cib = f"{simsuite_path}g_ocmb_cib5.fits" + + print(" ::: reading lensed cmb + cib + tsz090 map for hres:", h_cmb_cib_tsz090) + print(" ::: reading lensed cmb + cib + tsz150 map for hres:", h_cmb_cib_tsz150) + print(" ::: reading lensed cmb + cib map for grad:", g_cmb_cib) + + h_cmb_cib_tsz090_map = enmap.read_map(h_cmb_cib_tsz090, delayed=False) + h_cmb_cib_tsz150_map = enmap.read_map(h_cmb_cib_tsz150, delayed=False) + g_cmb_cib_map = enmap.read_map(g_cmb_cib, delayed=False) + +if args.cmb_ksz_msub: + + cmb_msub090 = f"{simsuite_path}h_ocmb_msub090_ksz_updated.fits" + cmb_msub150 = f"{simsuite_path}h_ocmb_msub150_ksz_updated.fits" + g_cmb_ksz = f"{simsuite_path}g_ocmb_ksz.fits" + + + print(" ::: reading the model subtracted maps!", cmb_msub090) + print(" ::: reading the model subtracted maps!", cmb_msub150) + print(" ::: reading lensed cmb + ksz map for grad:", g_cmb_ksz) + + h_cmb_ksz_msub090_map = enmap.read_map(cmb_msub090, delayed=False) + h_cmb_ksz_msub150_map = enmap.read_map(cmb_msub150, delayed=False) + g_cmb_ksz_map = enmap.read_map(g_cmb_ksz, delayed=False) + + + +if args.cmb_cib_msub: + + cmb_msub090 = f"{simsuite_path}h_ocmb_msub090_cib5_updated.fits" + cmb_msub150 = f"{simsuite_path}h_ocmb_msub150_cib5_updated.fits" + g_cmb_cib = f"{simsuite_path}g_ocmb_cib5.fits" + + print(" ::: reading the model subtracted maps!", cmb_msub090) + print(" ::: reading the model subtracted maps!", cmb_msub150) + print(" ::: reading lensed cmb + cib map for grad:", g_cmb_cib) + + h_cmb_cib_msub090_map = enmap.read_map(cmb_msub090, delayed=False) + h_cmb_cib_msub150_map = enmap.read_map(cmb_msub150, delayed=False) + g_cmb_cib_map = enmap.read_map(g_cmb_cib, delayed=False) + + + + +if args.cmb_ksz_cib: + + h_cmb_ksz_cib090 = f"{simsuite_path}h_ocmb_ksz_cib090_5.fits" + h_cmb_ksz_cib150 = f"{simsuite_path}h_ocmb_ksz_cib150_5.fits" + g_cmb_ksz_cib = f"{simsuite_path}g_ocmb_ksz_cib5.fits" + + print(" ::: reading lensed cmb + ksz + cib090 map for hres:", h_cmb_ksz_cib090) + print(" ::: reading lensed cmb + ksz + cib150 map for hres:", h_cmb_ksz_cib150) + print(" ::: reading lensed cmb + ksz + cib map for grad:", g_cmb_ksz_cib) + + h_cmb_ksz_cib090_map = enmap.read_map(h_cmb_ksz_cib090, delayed=False) + h_cmb_ksz_cib150_map = enmap.read_map(h_cmb_ksz_cib150, delayed=False) + g_cmb_ksz_cib_map = enmap.read_map(g_cmb_ksz_cib, delayed=False) + +if args.cmb_ksz_cib_tsz: + + h_cmb_tsz_ksz_cib090 = f"{simsuite_path}h_ocmb_tsz090_ksz_cib5.fits" + h_cmb_tsz_ksz_cib150 = f"{simsuite_path}h_ocmb_tsz150_ksz_cib5.fits" + g_cmb_ksz_cib = f"{simsuite_path}g_ocmb_ksz_cib5.fits" + + print(" ::: reading lensed cmb + tsz090 + ksz + cib090 map for hres:", h_cmb_tsz_ksz_cib090) + print(" ::: reading lensed cmb + tsz150 + ksz + cib150 map for hres:", h_cmb_tsz_ksz_cib150) + print(" ::: reading lensed cmb + ksz + cib map for grad:", g_cmb_ksz_cib) + + h_cmb_tsz_ksz_cib090_map = enmap.read_map(h_cmb_tsz_ksz_cib090, delayed=False) + h_cmb_tsz_ksz_cib150_map = enmap.read_map(h_cmb_tsz_ksz_cib150, delayed=False) + g_cmb_ksz_cib_map = enmap.read_map(g_cmb_ksz_cib, delayed=False) + +if args.cmb_ksz_cib_msub: + + cmb_msub090 = f"{simsuite_path}h_ocmb_msub090_ksz_cib5_updated.fits" + cmb_msub150 = f"{simsuite_path}h_ocmb_msub150_ksz_cib5_updated.fits" + g_cmb_ksz_cib = f"{simsuite_path}g_ocmb_ksz_cib5.fits" + + + print(" ::: reading the model subtracted maps!", cmb_msub090) + print(" ::: reading the model subtracted maps!", cmb_msub150) + print(" ::: reading lensed cmb + ksz + cib map for grad:", g_cmb_ksz_cib) + + h_cmb_ksz_cib090_msub_map = enmap.read_map(cmb_msub090, delayed=False) + h_cmb_ksz_cib150_msub_map = enmap.read_map(cmb_msub150, delayed=False) + g_cmb_ksz_cib_map = enmap.read_map(g_cmb_ksz_cib, delayed=False) - print(" ::: model image subtraction!") - cmb_msub_map090 = cmb_tsz_map090 - model_map090 - cmb_msub_map150 = cmb_tsz_map150 - model_map150 print(" ::: maps are ready!") @@ -353,8 +512,14 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): cov[:, 1, 0] = p12[sel] ms = np.stack([m1[sel], m2[sel]]).swapaxes(0, 1) rs = np.stack([b1[sel], b2[sel]]).swapaxes(0, 1) - num = np.linalg.solve(cov, ms) - den = np.linalg.solve(cov, rs) + + ## NumPy 2.x no longer supports implicit batched solve for (n,2,2)x(n,2), + ## so we explicitly solve each 2x2 system per ell. + # num = np.linalg.solve(cov, ms) + # den = np.linalg.solve(cov, rs) + num = np.array([np.linalg.solve(cov[i], ms[i]) for i in range(nells)]) + den = np.array([np.linalg.solve(cov[i], rs[i]) for i in range(nells)]) + tcov = 1.0 / np.einsum("ij,ij->i", rs, den) ksolve = np.einsum("ij,ij->i", rs, num) * tcov assert np.all(np.isfinite(ksolve)) @@ -383,6 +548,8 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): # this returns conc = 3.326 using Klypin et al 2016 c(M,z) # k2dmap is our fixed template on 2D map - this is tapered and filtered _,_,_,_,_,_,_,_,k2dmap,_ = lensing.kappa_nfw_profiley(mass=M200c,conc=None,z=z,z_s=1100.,background='critical',delta=200,apply_filter=True,lmin=lmin,lmax=lmax,res=res,rstamp=rstamp,rmin=rmin,rmax=rmax,rwidth=rwidth) +# # print(np.shape(k2dmap)) +# # io.plot_img(k2dmap, 'kplot.png') ##### TEST ##### --------------------------------------------------------------- @@ -404,7 +571,6 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): i = task cper = int((j + 1) / len(my_tasks) * 100.0) if rank==0: print(f"Rank {rank} performing task {task} as index {j} ({cper}% complete.).") - coords = np.array([decs[i], ras[i]]) * utils.degree @@ -440,7 +606,7 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): ) # lensed cmb for grad if args.cmb_tsz: - cmb_tsz150 = reproject.thumbnails( + h_cmb_tsz150 = reproject.thumbnails( cmb_tsz_map150, coords, r=maxr, @@ -450,7 +616,7 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): pixwin=False ) # lensed cmb + tsz150 - cmb_tsz090 = reproject.thumbnails( + h_cmb_tsz090 = reproject.thumbnails( cmb_tsz_map090, coords, r=maxr, @@ -461,7 +627,7 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): ) # lensed cmb + tsz090 if args.cmb_msub: - cmb_msub150 = reproject.thumbnails( + h_cmb_msub150 = reproject.thumbnails( cmb_msub_map150, coords, r=maxr, @@ -471,7 +637,7 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): pixwin=False ) # lensed cmb + tsz150 - model150 - cmb_msub090 = reproject.thumbnails( + h_cmb_msub090 = reproject.thumbnails( cmb_msub_map090, coords, r=maxr, @@ -481,6 +647,282 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): pixwin=False ) # lensed cmb + tsz090 - model090 + if args.cmb_ksz: + h_cmb_ksz = reproject.thumbnails( + h_cmb_ksz_map, + coords, + r=maxr, + res=px * utils.arcmin, + proj="tan", + oversample=2, + pixwin=False + ) # lensed cmb + ksz for hres + + g_cmb_ksz = reproject.thumbnails( + g_cmb_ksz_map, + coords, + r=maxr, + res=px * utils.arcmin, + proj="tan", + oversample=2, + pixwin=False + ) # lensed cmb + ksz for grad + + if args.cmb_cib: + h_cmb_cib090 = reproject.thumbnails( + h_cmb_cib090_map, + coords, + r=maxr, + res=px * utils.arcmin, + proj="tan", + oversample=2, + pixwin=False + ) # lensed cmb + cib090 for hres + + h_cmb_cib150 = reproject.thumbnails( + h_cmb_cib150_map, + coords, + r=maxr, + res=px * utils.arcmin, + proj="tan", + oversample=2, + pixwin=False + ) # lensed cmb + cib150 for hres + + g_cmb_cib = reproject.thumbnails( + g_cmb_cib_map, + coords, + r=maxr, + res=px * utils.arcmin, + proj="tan", + oversample=2, + pixwin=False + ) # lensed cmb + cib150 for grad + + + + + if args.cmb_ksz_tsz: + h_cmb_ksz_tsz090 = reproject.thumbnails( + h_cmb_ksz_tsz090_map, + coords, + r=maxr, + res=px * utils.arcmin, + proj="tan", + oversample=2, + pixwin=False + ) # lensed cmb + ksz + tsz090 for hres + + h_cmb_ksz_tsz150 = reproject.thumbnails( + h_cmb_ksz_tsz150_map, + coords, + r=maxr, + res=px * utils.arcmin, + proj="tan", + oversample=2, + pixwin=False + ) # lensed cmb + ksz + tsz150 for hres + + g_cmb_ksz = reproject.thumbnails( + g_cmb_ksz_map, + coords, + r=maxr, + res=px * utils.arcmin, + proj="tan", + oversample=2, + pixwin=False + ) # lensed cmb + ksz for grad + + if args.cmb_cib_tsz: + h_cmb_cib_tsz090 = reproject.thumbnails( + h_cmb_cib_tsz090_map, + coords, + r=maxr, + res=px * utils.arcmin, + proj="tan", + oversample=2, + pixwin=False + ) # lensed cmb + cib + tsz090 for hres + + h_cmb_cib_tsz150 = reproject.thumbnails( + h_cmb_cib_tsz150_map, + coords, + r=maxr, + res=px * utils.arcmin, + proj="tan", + oversample=2, + pixwin=False + ) # lensed cmb + cib + tsz150 for hres + + g_cmb_cib = reproject.thumbnails( + g_cmb_cib_map, + coords, + r=maxr, + res=px * utils.arcmin, + proj="tan", + oversample=2, + pixwin=False + ) # lensed cmb + cib for grad + + if args.cmb_ksz_msub: + h_cmb_ksz_msub090 = reproject.thumbnails( + h_cmb_ksz_msub090_map, + coords, + r=maxr, + res=px * utils.arcmin, + proj="tan", + oversample=2, + pixwin=False + ) # lensed cmb + ksz + tsz090 - model for hres + + h_cmb_ksz_msub150 = reproject.thumbnails( + h_cmb_ksz_msub150_map, + coords, + r=maxr, + res=px * utils.arcmin, + proj="tan", + oversample=2, + pixwin=False + ) # lensed cmb + ksz + tsz150 - model for hres + + g_cmb_ksz = reproject.thumbnails( + g_cmb_ksz_map, + coords, + r=maxr, + res=px * utils.arcmin, + proj="tan", + oversample=2, + pixwin=False + ) # lensed cmb + ksz for grad + + if args.cmb_cib_msub: + h_cmb_cib_msub090 = reproject.thumbnails( + h_cmb_cib_msub090_map, + coords, + r=maxr, + res=px * utils.arcmin, + proj="tan", + oversample=2, + pixwin=False + ) # lensed cmb + cib + tsz090 - model for hres + + h_cmb_cib_msub150 = reproject.thumbnails( + h_cmb_cib_msub150_map, + coords, + r=maxr, + res=px * utils.arcmin, + proj="tan", + oversample=2, + pixwin=False + ) # lensed cmb + cib + tsz150 - model for hres + + g_cmb_cib = reproject.thumbnails( + g_cmb_cib_map, + coords, + r=maxr, + res=px * utils.arcmin, + proj="tan", + oversample=2, + pixwin=False + ) # lensed cmb + cib for grad + + + + + + + if args.cmb_ksz_cib: + h_cmb_ksz_cib090 = reproject.thumbnails( + h_cmb_ksz_cib090_map, + coords, + r=maxr, + res=px * utils.arcmin, + proj="tan", + oversample=2, + pixwin=False + ) # lensed cmb + ksz + cib090 for hres + + h_cmb_ksz_cib150 = reproject.thumbnails( + h_cmb_ksz_cib150_map, + coords, + r=maxr, + res=px * utils.arcmin, + proj="tan", + oversample=2, + pixwin=False + ) # lensed cmb + ksz + cib150 for hres + + g_cmb_ksz_cib = reproject.thumbnails( + g_cmb_ksz_cib_map, + coords, + r=maxr, + res=px * utils.arcmin, + proj="tan", + oversample=2, + pixwin=False + ) # lensed cmb + ksz + cib150 for grad + + if args.cmb_ksz_cib_tsz: + h_cmb_tsz_ksz_cib090 = reproject.thumbnails( + h_cmb_tsz_ksz_cib090_map, + coords, + r=maxr, + res=px * utils.arcmin, + proj="tan", + oversample=2, + pixwin=False + ) # lensed cmb + tsz090 + ksz + cib090 for hres + + h_cmb_tsz_ksz_cib150 = reproject.thumbnails( + h_cmb_tsz_ksz_cib150_map, + coords, + r=maxr, + res=px * utils.arcmin, + proj="tan", + oversample=2, + pixwin=False + ) # lensed cmb + tsz150 + ksz + cib150 for hres + + g_cmb_ksz_cib = reproject.thumbnails( + g_cmb_ksz_cib_map, + coords, + r=maxr, + res=px * utils.arcmin, + proj="tan", + oversample=2, + pixwin=False + ) # lensed cmb + ksz + cib for grad + + if args.cmb_ksz_cib_msub: + h_cmb_ksz_cib090_msub = reproject.thumbnails( + h_cmb_ksz_cib090_msub_map, + coords, + r=maxr, + res=px * utils.arcmin, + proj="tan", + oversample=2, + pixwin=False + ) # lensed cmb + tsz090 + ksz + cib090 - model for hres + + h_cmb_ksz_cib150_msub = reproject.thumbnails( + h_cmb_ksz_cib150_msub_map, + coords, + r=maxr, + res=px * utils.arcmin, + proj="tan", + oversample=2, + pixwin=False + ) # lensed cmb + tsz150 + ksz + cib150 - model for hres + + g_cmb_ksz_cib = reproject.thumbnails( + g_cmb_ksz_cib_map, + coords, + r=maxr, + res=px * utils.arcmin, + proj="tan", + oversample=2, + pixwin=False + ) # lensed cmb + ksz + cib for grad # initialise calculations based on geometry if j == 0: @@ -495,12 +937,57 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): assert wcsutils.equal(kstamp.wcs, g_cmb.wcs) if args.cmb_tsz: - assert wcsutils.equal(kstamp.wcs, cmb_tsz150.wcs) - assert wcsutils.equal(kstamp.wcs, cmb_tsz090.wcs) - - if args.cmb_msub: - assert wcsutils.equal(kstamp.wcs, cmb_msub150.wcs) - assert wcsutils.equal(kstamp.wcs, cmb_msub090.wcs) + assert wcsutils.equal(kstamp.wcs, h_cmb_tsz150.wcs) + assert wcsutils.equal(kstamp.wcs, h_cmb_tsz090.wcs) + + elif args.cmb_msub: + assert wcsutils.equal(kstamp.wcs, h_cmb_msub150.wcs) + assert wcsutils.equal(kstamp.wcs, h_cmb_msub090.wcs) + + elif args.cmb_ksz: + assert wcsutils.equal(kstamp.wcs, h_cmb_ksz.wcs) + assert wcsutils.equal(kstamp.wcs, g_cmb_ksz.wcs) + + elif args.cmb_cib: + assert wcsutils.equal(kstamp.wcs, h_cmb_cib150.wcs) + assert wcsutils.equal(kstamp.wcs, h_cmb_cib090.wcs) + assert wcsutils.equal(kstamp.wcs, g_cmb_cib.wcs) + + elif args.cmb_ksz_tsz: + assert wcsutils.equal(kstamp.wcs, h_cmb_ksz_tsz150.wcs) + assert wcsutils.equal(kstamp.wcs, h_cmb_ksz_tsz090.wcs) + assert wcsutils.equal(kstamp.wcs, g_cmb_ksz.wcs) + + elif args.cmb_cib_tsz: + assert wcsutils.equal(kstamp.wcs, h_cmb_cib_tsz150.wcs) + assert wcsutils.equal(kstamp.wcs, h_cmb_cib_tsz090.wcs) + assert wcsutils.equal(kstamp.wcs, g_cmb_cib.wcs) + + elif args.cmb_ksz_msub: + assert wcsutils.equal(kstamp.wcs, h_cmb_ksz_msub150.wcs) + assert wcsutils.equal(kstamp.wcs, h_cmb_ksz_msub090.wcs) + assert wcsutils.equal(kstamp.wcs, g_cmb_ksz.wcs) + + elif args.cmb_cib_msub: + assert wcsutils.equal(kstamp.wcs, h_cmb_cib_msub150.wcs) + assert wcsutils.equal(kstamp.wcs, h_cmb_cib_msub090.wcs) + assert wcsutils.equal(kstamp.wcs, g_cmb_cib.wcs) + + + elif args.cmb_ksz_cib: + assert wcsutils.equal(kstamp.wcs, h_cmb_ksz_cib150.wcs) + assert wcsutils.equal(kstamp.wcs, h_cmb_ksz_cib090.wcs) + assert wcsutils.equal(kstamp.wcs, g_cmb_ksz_cib.wcs) + + elif args.cmb_tsz: + assert wcsutils.equal(kstamp.wcs, h_cmb_tsz_ksz_cib150.wcs) + assert wcsutils.equal(kstamp.wcs, h_cmb_tsz_ksz_cib090.wcs) + assert wcsutils.equal(kstamp.wcs, g_cmb_ksz_cib.wcs) + + elif args.cmb_ksz_cib_msub: + assert wcsutils.equal(kstamp.wcs, h_cmb_ksz_cib150_msub.wcs) + assert wcsutils.equal(kstamp.wcs, h_cmb_ksz_cib090_msub.wcs) + assert wcsutils.equal(kstamp.wcs, g_cmb_ksz_cib.wcs) # get an edge taper map and apodize taper = maps.get_taper( @@ -516,7 +1003,9 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): beam2d = maps.gauss_beam(modlmap, fwhm) beam2d_plc = maps.gauss_beam(modlmap, fwhm_plc) - if args.cmb_tsz or args.cmb_msub: + # if args.cmb_tsz or args.cmb_msub or args.cmb_cib or args.cmb_ksz_cib or args.cmb_ksz_cib_tsz or args.cmb_ksz_cib_msub: + if not (args.cmb or args.cmb_ksz): + bfunc150 = lambda x: maps.gauss_beam(ilc_beam_fwhm150, x) bfunc90 = lambda x: maps.gauss_beam(ilc_beam_fwhm090, x) @@ -533,7 +1022,7 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): theory = cosmology.default_theory() ucltt2d = theory.lCl("TT", modlmap) - if args.cmb_tsz or args.cmb_msub: + if not (args.cmb or args.cmb_ksz): # bin size and range for 1D binned power spectrum minell = 2 * maps.minimum_ell(shape, wcs) l_edges = np.arange(minell / 2, 8001, minell) @@ -543,39 +1032,86 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): # fsky for bandpower variance fsky = enmap.area(shape, wcs) * w2 / 4.0 / np.pi - # same filter as the post-reconstuction for true kappa (also tapered!) k_stamp = maps.filter_map(kstamp*taper, kmask) s.add_to_stack("kstamp", k_stamp) binned_true = bin(k_stamp, modrmap * (180 * 60 / np.pi), bin_edges) s.add_to_stats("tk1d", binned_true) - - if args.cmb_tsz: - hres150 = cmb_tsz150 - hres90 = cmb_tsz090 - elif args.cmb_msub: - hres150 = cmb_msub150 - hres90 = cmb_msub090 - else: + if args.cmb: hres = h_cmb - - - # to filter weird stamps when model image subtraction is done for the mean-field - if args.cmb_msub: + g_cmb = g_cmb + elif args.cmb_tsz: + hres150 = h_cmb_tsz150 + hres90 = h_cmb_tsz090 + g_cmb = g_cmb + elif args.cmb_msub: + hres150 = h_cmb_msub150 + hres90 = h_cmb_msub090 + g_cmb = g_cmb + elif args.cmb_ksz: + hres = h_cmb_ksz + g_cmb = g_cmb_ksz + elif args.cmb_cib: + hres150 = h_cmb_cib150 + hres90 = h_cmb_cib090 + g_cmb = g_cmb_cib + + elif args.cmb_ksz_tsz: + hres150 = h_cmb_ksz_tsz150 + hres90 = h_cmb_ksz_tsz090 + g_cmb = g_cmb_ksz + elif args.cmb_cib_tsz: + hres150 = h_cmb_cib_tsz150 + hres90 = h_cmb_cib_tsz090 + g_cmb = g_cmb_cib + elif args.cmb_ksz_msub: + hres150 = h_cmb_ksz_msub150 + hres90 = h_cmb_ksz_msub090 + g_cmb = g_cmb_ksz + elif args.cmb_cib_msub: + hres150 = h_cmb_cib_msub150 + hres90 = h_cmb_cib_msub090 + g_cmb = g_cmb_cib + + + elif args.cmb_ksz_cib: + hres150 = h_cmb_ksz_cib150 + hres90 = h_cmb_ksz_cib090 + g_cmb = g_cmb_ksz_cib + elif args.cmb_ksz_cib_tsz: + hres150 = h_cmb_tsz_ksz_cib150 + hres90 = h_cmb_tsz_ksz_cib090 + g_cmb = g_cmb_ksz_cib + elif args.cmb_ksz_cib_msub: + hres150 = h_cmb_ksz_cib150_msub + hres90 = h_cmb_ksz_cib090_msub + g_cmb = g_cmb_ksz_cib + + + # filter weird stamps (temporary1) + if not (args.cmb or args.cmb_ksz): if np.any(abs(hres150) > 1e3) or np.any(abs(hres90) > 1e3): print(f"{task} has anomalously high hres150 or hres90") continue + else: + if np.any(abs(hres) > 1e3): + print(f"{task} has anomalously high hres") + continue # taper stamp tapered_grad = g_cmb * taper - if args.cmb_tsz or args.cmb_msub: + if not (args.cmb or args.cmb_ksz): tapered_hres150 = hres150 * taper tapered_hres90 = hres90 * taper else: tapered_hres = hres * taper + + + + s.add_to_stack("grad2d_before", g_cmb) g_filtered = maps.filter_map(g_cmb, ymask) s.add_to_stack("grad2d_before_filtered", g_filtered) @@ -583,6 +1119,16 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): s.add_to_stats("grad1d_before", binned_grad) + # # this is a new inpainting method + # if args.inpaint: + # mask = np.zeros(g_cmb.shape, dtype=bool) + # mask[g_cmb.modrmap()<10.0 * utils.arcmin] = True + # mask[mask==0] = False + # mask = enmap.enmap(mask, g_cmb.wcs) + # g_cmb = maps.gapfill_edge_conv_flat(g_cmb, mask) + # tapered_grad = g_cmb * taper + + # this is the old inpainting methood if args.inpaint: # for gradient leg (should be used with args.is_observed at the moment) """ @@ -607,7 +1153,7 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): pcov = pixcov.tpcov_from_ivar(Ndown, pivar, theory.lCl, beam_fn) geo = pixcov.make_geometry(pshape, pwcs, rmin, n=Ndown, deproject=True, iau=False, res=None, pcov=pcov) - cutout = pixcov.inpaint_stamp(cutout, geo) + cutout = pixcov.inpaint_stamp(cutout, geo, add_noise=False) tapered_grad[cutout_sel] = cutout.copy() tapered_grad = np.nan_to_num(tapered_grad) g_cmb = tapered_grad / taper @@ -622,8 +1168,7 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): - if args.cmb_tsz or args.cmb_msub: - + if not (args.cmb or args.cmb_ksz): # ILC will return a beam-deconvolved Fourier transformed stamp k_grad = enmap.fft(tapered_grad, normalize="phys") k150 = enmap.fft(tapered_hres150, normalize="phys") @@ -759,7 +1304,7 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): assert np.all(np.isfinite(feed_dict[key])) # ask for reconstruction in Fourier space - cqe = symlens.QE( + cqe = QE( shape, wcs, feed_dict, @@ -831,19 +1376,56 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): # transform to real space kappa = enmap.ifft(rkmap, normalize="phys").real + # filter weird stamps (temporary2) + if np.any(np.abs(kappa) > 15): + print(f"{task} has large kappa") + continue + + # stack reconstructed kappa s.add_to_stack("lstamp", kappa) binned_kappa = bin(kappa, modrmap * (180 * 60 / np.pi), bin_edges) s.add_to_stats("k1d", binned_kappa) # stack stamps pre-reconstruction as well (same filter as hres leg) + if args.cmb: + hres = h_cmb + grad = g_cmb if args.cmb_tsz: - hres = cmb_tsz150 + hres = h_cmb_tsz150 + grad = g_cmb elif args.cmb_msub: - hres = cmb_msub150 - else: - hres = h_cmb - grad = g_cmb + hres = h_cmb_msub150 + grad = g_cmb + elif args.cmb_ksz: + hres = h_cmb_ksz + grad = g_cmb_ksz + elif args.cmb_cib: + hres = h_cmb_cib150 + grad = g_cmb_cib + + elif args.cmb_ksz_tsz: + hres = h_cmb_ksz_tsz150 + grad = g_cmb_ksz + elif args.cmb_cib_tsz: + hres = h_cmb_cib_tsz150 + grad = g_cmb_cib + elif args.cmb_ksz_msub: + hres = h_cmb_ksz_msub150 + grad = g_cmb_ksz + elif args.cmb_cib_msub: + hres = h_cmb_cib_msub150 + grad = g_cmb_cib + + elif args.cmb_ksz_cib: + hres = h_cmb_ksz_cib150 + grad = g_cmb_ksz_cib + elif args.cmb_ksz_cib_tsz: + hres = h_cmb_tsz_ksz_cib150 + grad = g_cmb_ksz_cib + elif args.cmb_ksz_cib_msub: + hres = h_cmb_ksz_cib150_msub + grad = g_cmb_ksz_cib hres_stamp = maps.filter_map(hres, ymask) grad_stamp = maps.filter_map(grad, ymask) diff --git a/stack.py b/stack.py index a60c781..10ad1b4 100644 --- a/stack.py +++ b/stack.py @@ -1,6 +1,3 @@ -# inpainting test -# do inpainting on 90 and 150 GHz separately first and then ILC - import numpy as np import matplotlib # Force matplotlib to not use any Xwindows backend. @@ -170,7 +167,7 @@ rms_map = maps.rms_from_ivar( - imap_90, cylindrical=True + imap_90, cylindrical=True, safe=False # safe=False for DR6 maps ) # convert to RMS noise map @@ -682,42 +679,56 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): else: gact_stamp_150 = gastamp_150 * taper gact_stamp_90 = gastamp_90 * taper - - if args.inpaint: # only with hres_grad for now - """ - If inpainting, we - (1) resample the stamp to 64x64 (2 arcmin pixels) - (2) Inpaint a hole of radius 4 arcmin - """ - rmin = 4 * utils.arcmin - crop_pixels = int(16. / args.pwidth) # 16 arcminutes wide - act150 = maps.crop_center(gact_stamp_150,cropy=crop_pixels,cropx=crop_pixels,sel=False) - act90 = maps.crop_center(gact_stamp_90,cropy=crop_pixels,cropx=crop_pixels,sel=False) - act_sel150 = maps.crop_center(gact_stamp_150,cropy=crop_pixels,cropx=crop_pixels,sel=True) - act_sel90 = maps.crop_center(gact_stamp_90,cropy=crop_pixels,cropx=crop_pixels,sel=True) - Ndown150, Ndown2 = act150.shape[-2:] - #print(Ndown150, Ndown2) # 32 32 - if Ndown150 != Ndown2: raise Exception - Ndown90, Ndown2 = act90.shape[-2:] - if Ndown90 != Ndown2: raise Exception + + + + + + if args.inpaint: # here inpainting is done on gradient leg map + + s.add_to_stack("grad2d_before", pstamp) + binned_grad = bin(pstamp, pstamp.modrmap() * (180 * 60 / np.pi), bin_edges) + s.add_to_stats("grad1d_before", binned_grad) + + # # this is a NEW inpainting method + # mask = np.zeros(pstamp.shape, dtype=bool) + # mask[pstamp.modrmap()<10.0 * utils.arcmin] = True + # mask[mask==0] = False + # mask = enmap.enmap(mask, pstamp.wcs) + # pstamp = maps.gapfill_edge_conv_flat(pstamp, mask) + # plc_stamp = pstamp * taper + + # this is a OLD inpainting method + rmin = 4.0 * utils.arcmin # inpainting radius is set to 4 + px = 0.5 + crop_pixels = int(16. / px) # 16 arcminutes wide + tapered_grad = plc_stamp + cutout = maps.crop_center(tapered_grad, cropy=crop_pixels, cropx=crop_pixels, sel=False) + cutout_sel = maps.crop_center(tapered_grad, cropy=crop_pixels, cropx=crop_pixels, sel=True) + Ndown, Ndown2 = cutout.shape[-2:] + if Ndown != Ndown2: raise Exception + grad_fiducial_rms = 35 if j==0: from orphics import pixcov - pshape = act150.shape - pwcs = act150.wcs - beam_fn150 = cutils.load_beam("f150") - beam_fn90 = cutils.load_beam("f090") + pshape = cutout.shape + pwcs = cutout.wcs + fwhm = 5.0 + beam_fn = lambda x: maps.gauss_beam(fwhm, x) ipsizemap = enmap.pixsizemap(pshape, pwcs) - pivar = maps.ivar(pshape, pwcs, defaults.highres_fiducial_rms, ipsizemap=ipsizemap) - pcov150 = pixcov.tpcov_from_ivar(Ndown150, pivar, theory.lCl, beam_fn150) - pcov90 = pixcov.tpcov_from_ivar(Ndown90, pivar, theory.lCl, beam_fn90) - geo150 = pixcov.make_geometry(pshape, pwcs, rmin, n=Ndown150, deproject=True, iau=False, res=None, pcov=pcov150) - geo90 = pixcov.make_geometry(pshape, pwcs, rmin, n=Ndown90, deproject=True, iau=False, res=None, pcov=pcov90) - - act150 = pixcov.inpaint_stamp(act150, geo150) - act90 = pixcov.inpaint_stamp(act90, geo90) - gact_stamp_150[act_sel150] = act150.copy() - gact_stamp_90[act_sel90] = act90.copy() + pivar = maps.ivar(pshape, pwcs, grad_fiducial_rms, ipsizemap=ipsizemap) + pcov = pixcov.tpcov_from_ivar(Ndown, pivar, theory.lCl, beam_fn) + geo = pixcov.make_geometry(pshape, pwcs, rmin, n=Ndown, deproject=True, iau=False, res=None, pcov=pcov) + + cutout = pixcov.inpaint_stamp(cutout, geo) + tapered_grad[cutout_sel] = cutout.copy() + plc_stamp = tapered_grad + pstamp = tapered_grad / taper + + s.add_to_stack("grad2d_after", pstamp) + binned_grad = bin(pstamp, pstamp.modrmap() * (180 * 60 / np.pi), bin_edges) + s.add_to_stats("grad1d_after", binned_grad) + if args.day_null: nact_stamp_150 = nastamp_150 * taper From f1a2cdbaba25eab715edab492fa15ed013112a6d Mon Sep 17 00:00:00 2001 From: Mathew Madhavacheril Date: Tue, 19 May 2026 19:30:54 -0400 Subject: [PATCH 05/14] more robust inpainting --- inpaint_utils.py | 152 +++++++++++++++++++++++++++++++++++++++++++++++ recon_sim.py | 55 +++++------------ stack.py | 73 ++++++++--------------- 3 files changed, 191 insertions(+), 89 deletions(-) create mode 100644 inpaint_utils.py diff --git a/inpaint_utils.py b/inpaint_utils.py new file mode 100644 index 0000000..b150d6c --- /dev/null +++ b/inpaint_utils.py @@ -0,0 +1,152 @@ +"""Gradient-leg inpainting helpers for cmbhalolensing. + +Provides `GradientInpainter`, a class that encapsulates pixel-pixel +covariance inpainting on a coarse sub-stamp (60'x60' / 2' by default), +then upgrades and substitutes only the hole pixels back into the +fine-resolution input. Built once per channel at the first cluster, +called per cluster after that. + +Designed to be shared between `recon_sim.py` (sims) and `stack.py` +(data). For data with per-frequency inpainting before ILC, instantiate +one `GradientInpainter` per frequency with different `beam_fn`. +""" + +import numpy as np +from pixell import enmap, utils +from orphics import maps, pixcov + + +class Inpainter: + """Pixcov-based inpainting on a coarse sub-stamp. + + One instance binds a fixed (stamp geometry, beam, white noise, + hole radius). + + Parameters + ---------- + fine_shape, fine_wcs + Geometry of the fine-resolution stamp that will be passed to + ``inpaint()``. Used to build the fine-grid hole mask and the + cutout selector. + hole_radius_arcmin + Radius of the circular hole. + theory + Object providing ``theory.lCl('TT', ell)``. + beam_fn + Callable ``beam_fn(ell) -> B(ell)``. For Gaussian beams, + ``lambda x: maps.gauss_beam(fwhm, x)``. For per-frequency + tabulated beams (data), pass the loader's callable directly. + noise_uK_arcmin + White-noise level for the diagonal noise covariance. + context_width_arcmin + Total width of the coarse sub-stamp. Default 60'. Must + comfortably exceed 2 * hole_radius_arcmin. + pcov_res_arcmin + Resolution of the coarse pcov grid. Default 2'. + ``context_width_arcmin / pcov_res_arcmin`` must be an integer + multiple of the fine pixel size. + eigval_floor + Floor on the conditional-covariance eigenvalues, relative to + the largest. Only used when ``inpaint(add_noise=True)``; + harmless otherwise. + """ + + def __init__( + self, + fine_shape, + fine_wcs, + hole_radius_arcmin, + theory, + beam_fn, + noise_uK_arcmin, + context_width_arcmin=60.0, + pcov_res_arcmin=2.0, + eigval_floor=1e-6, + ): + fine_res = np.abs(fine_wcs.wcs.cdelt[1]) * 60.0 # arcmin + self._dfac = int(round(pcov_res_arcmin / fine_res)) + self._fine_crop = int(round(context_width_arcmin / fine_res)) + coarse_n = self._fine_crop // self._dfac + if coarse_n * self._dfac != self._fine_crop: + raise ValueError( + "context_width_arcmin must be an integer multiple of " + "pcov_res_arcmin AND of fine_res " + f"(fine_res={fine_res}', context={context_width_arcmin}', " + f"pcov_res={pcov_res_arcmin}')" + ) + if 2 * hole_radius_arcmin >= context_width_arcmin: + raise ValueError( + f"context_width_arcmin ({context_width_arcmin}) must exceed " + f"2 * hole_radius_arcmin ({2 * hole_radius_arcmin})" + ) + + # Coarse pcov geometry built from a dummy fine map so its WCS + # stays in the thumbnail's TAN projection. + dummy = enmap.zeros(fine_shape, fine_wcs) + cut0 = maps.crop_center( + dummy, cropy=self._fine_crop, cropx=self._fine_crop, sel=False + ) + coarse0 = enmap.downgrade(cut0, self._dfac) + pshape, pwcs = coarse0.shape, coarse0.wcs + + # Diagonal white-noise + theory + beam -> pixcov on coarse grid + ipsizemap = enmap.pixsizemap(pshape, pwcs) + pivar = maps.ivar(pshape, pwcs, noise_uK_arcmin, ipsizemap=ipsizemap) + pcov = pixcov.tpcov_from_ivar(coarse_n, pivar, theory.lCl, beam_fn) + + hole_radius = hole_radius_arcmin * utils.arcmin + self._geo = pixcov.make_geometry( + pshape, + pwcs, + hole_radius, + n=coarse_n, + deproject=True, + iau=False, + res=None, + pcov=pcov, + eigval_floor=eigval_floor, + ) + + self._cutout_sel = maps.crop_center( + dummy, cropy=self._fine_crop, cropx=self._fine_crop, sel=True + ) + modrmap = enmap.modrmap(fine_shape, fine_wcs) + self._fine_hole_mask = modrmap < hole_radius + self._box_hole_mask = self._fine_hole_mask[self._cutout_sel] + self._fine_shape = tuple(fine_shape) + + def inpaint(self, fine_stamp, add_noise=False, rng=None): + """Return a copy of ``fine_stamp`` with the central hole filled. + + Only the hole pixels are modified; surrounding pixels are left + exactly as in the input (no downgrade/upgrade smearing). + + Parameters + ---------- + fine_stamp : enmap + Fine-resolution stamp. Must match the (shape, wcs) the + inpainter was built with. + add_noise : bool + If True, add a constrained-realisation random draw on top of + the conditional mean. Default False = mean-only. + rng : np.random.Generator, optional + RNG for the random draw. Ignored when ``add_noise=False``. + """ + if tuple(fine_stamp.shape) != self._fine_shape: + raise ValueError( + f"fine_stamp shape {tuple(fine_stamp.shape)} != " + f"geometry shape {self._fine_shape}" + ) + + out = fine_stamp.copy() + fine_cutout = maps.crop_center( + out, cropy=self._fine_crop, cropx=self._fine_crop, sel=False + ).copy() + coarse = enmap.downgrade(fine_cutout, self._dfac) + coarse_inp = pixcov.inpaint_stamp( + coarse, self._geo, add_noise=add_noise, rng=rng + ) + fine_inp = enmap.upgrade(coarse_inp, self._dfac) + fine_cutout[self._box_hole_mask] = fine_inp[self._box_hole_mask] + out[self._cutout_sel] = fine_cutout + return out diff --git a/recon_sim.py b/recon_sim.py index 2133f3f..f7dca0d 100644 --- a/recon_sim.py +++ b/recon_sim.py @@ -1099,8 +1099,6 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): print(f"{task} has anomalously high hres") continue - # taper stamp - tapered_grad = g_cmb * taper if not (args.cmb or args.cmb_ksz): tapered_hres150 = hres150 * taper @@ -1118,46 +1116,21 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): binned_grad = bin(g_cmb, g_cmb.modrmap() * (180 * 60 / np.pi), bin_edges) s.add_to_stats("grad1d_before", binned_grad) + if args.inpaint: + if j == 0: + from inpaint_utils import Inpainter + grad_inpainter = Inpainter( + fine_shape=g_cmb.shape, + fine_wcs=g_cmb.wcs, + hole_radius_arcmin=inpaint_radius, + theory=theory, + beam_fn=lambda x: maps.gauss_beam(fwhm_plc, x), + noise_uK_arcmin=nlevel_plc, + ) + g_cmb = grad_inpainter.inpaint(g_cmb) - # # this is a new inpainting method - # if args.inpaint: - # mask = np.zeros(g_cmb.shape, dtype=bool) - # mask[g_cmb.modrmap()<10.0 * utils.arcmin] = True - # mask[mask==0] = False - # mask = enmap.enmap(mask, g_cmb.wcs) - # g_cmb = maps.gapfill_edge_conv_flat(g_cmb, mask) - # tapered_grad = g_cmb * taper - - - # this is the old inpainting methood - if args.inpaint: # for gradient leg (should be used with args.is_observed at the moment) - """ - If inpainting, we - (1) resample the stamp to 64x64 (2 arcmin pixels) - (2) Inpaint a hole of radius 4 arcmin - """ - rmin = inpaint_radius * utils.arcmin - crop_pixels = int(inpaint_radius * 4. / px) # 16 arcminutes wide for 4 arcmin radius - cutout = maps.crop_center(tapered_grad, cropy=crop_pixels, cropx=crop_pixels, sel=False) - cutout_sel = maps.crop_center(tapered_grad, cropy=crop_pixels, cropx=crop_pixels, sel=True) - Ndown, Ndown2 = cutout.shape[-2:] - if Ndown != Ndown2: raise Exception - - if j==0: # set up for planck-like beam and noise - from orphics import pixcov - pshape = cutout.shape - pwcs = cutout.wcs - beam_fn = lambda x: maps.gauss_beam(fwhm_plc, x) - ipsizemap = enmap.pixsizemap(pshape, pwcs) - pivar = maps.ivar(pshape, pwcs, nlevel_plc, ipsizemap=ipsizemap) - pcov = pixcov.tpcov_from_ivar(Ndown, pivar, theory.lCl, beam_fn) - geo = pixcov.make_geometry(pshape, pwcs, rmin, n=Ndown, deproject=True, iau=False, res=None, pcov=pcov) - - cutout = pixcov.inpaint_stamp(cutout, geo, add_noise=False) - tapered_grad[cutout_sel] = cutout.copy() - tapered_grad = np.nan_to_num(tapered_grad) - g_cmb = tapered_grad / taper - g_cmb = np.nan_to_num(g_cmb) + # taper stamp + tapered_grad = g_cmb * taper s.add_to_stack("grad2d_after", g_cmb) diff --git a/stack.py b/stack.py index 10ad1b4..f8c9e4f 100644 --- a/stack.py +++ b/stack.py @@ -674,60 +674,37 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): act_stamp_150 = astamp_150 * taper act_stamp_90 = astamp_90 * taper + + if args.inpaint: # gradient leg (Planck-like, single channel) + if args.hres_grad: raise NotImplementedError("Inpainting when using high-resolution gradient not supported.") + + s.add_to_stack("grad2d_before", pstamp) + binned_grad = bin(pstamp, pstamp.modrmap() * (180 * 60 / np.pi), bin_edges) + s.add_to_stats("grad1d_before", binned_grad) + + if j == 0: + from inpaint_utils import Inpainter + grad_inpainter = Inpainter( + fine_shape=pstamp.shape, + fine_wcs=pstamp.wcs, + hole_radius_arcmin=4.0, + theory=theory, + beam_fn=lambda x: maps.gauss_beam(5.0, x), + noise_uK_arcmin=35.0, + ) + pstamp = grad_inpainter.inpaint(pstamp) + + s.add_to_stack("grad2d_after", pstamp) + binned_grad = bin(pstamp, pstamp.modrmap() * (180 * 60 / np.pi), bin_edges) + s.add_to_stats("grad1d_after", binned_grad) + if not (args.hres_grad): plc_stamp = pstamp * taper else: gact_stamp_150 = gastamp_150 * taper gact_stamp_90 = gastamp_90 * taper + - - - - - if args.inpaint: # here inpainting is done on gradient leg map - - s.add_to_stack("grad2d_before", pstamp) - binned_grad = bin(pstamp, pstamp.modrmap() * (180 * 60 / np.pi), bin_edges) - s.add_to_stats("grad1d_before", binned_grad) - - # # this is a NEW inpainting method - # mask = np.zeros(pstamp.shape, dtype=bool) - # mask[pstamp.modrmap()<10.0 * utils.arcmin] = True - # mask[mask==0] = False - # mask = enmap.enmap(mask, pstamp.wcs) - # pstamp = maps.gapfill_edge_conv_flat(pstamp, mask) - # plc_stamp = pstamp * taper - - # this is a OLD inpainting method - rmin = 4.0 * utils.arcmin # inpainting radius is set to 4 - px = 0.5 - crop_pixels = int(16. / px) # 16 arcminutes wide - tapered_grad = plc_stamp - cutout = maps.crop_center(tapered_grad, cropy=crop_pixels, cropx=crop_pixels, sel=False) - cutout_sel = maps.crop_center(tapered_grad, cropy=crop_pixels, cropx=crop_pixels, sel=True) - Ndown, Ndown2 = cutout.shape[-2:] - if Ndown != Ndown2: raise Exception - grad_fiducial_rms = 35 - - if j==0: - from orphics import pixcov - pshape = cutout.shape - pwcs = cutout.wcs - fwhm = 5.0 - beam_fn = lambda x: maps.gauss_beam(fwhm, x) - ipsizemap = enmap.pixsizemap(pshape, pwcs) - pivar = maps.ivar(pshape, pwcs, grad_fiducial_rms, ipsizemap=ipsizemap) - pcov = pixcov.tpcov_from_ivar(Ndown, pivar, theory.lCl, beam_fn) - geo = pixcov.make_geometry(pshape, pwcs, rmin, n=Ndown, deproject=True, iau=False, res=None, pcov=pcov) - - cutout = pixcov.inpaint_stamp(cutout, geo) - tapered_grad[cutout_sel] = cutout.copy() - plc_stamp = tapered_grad - pstamp = tapered_grad / taper - - s.add_to_stack("grad2d_after", pstamp) - binned_grad = bin(pstamp, pstamp.modrmap() * (180 * 60 / np.pi), bin_edges) - s.add_to_stats("grad1d_after", binned_grad) if args.day_null: From 36c2c1be12ca6506d63df2bf97c1991ef1c1a055 Mon Sep 17 00:00:00 2001 From: Mathew Madhavacheril Date: Tue, 19 May 2026 21:27:32 -0400 Subject: [PATCH 06/14] calibration --- cal_instructions.md | 87 ++++++++ cal_response.py | 434 ++++++++++++++++++++++++++++++++++++++++ cal_utils.py | 98 +++++++++ input/cal_recon_sim.yml | 33 +++ input/cal_stack.yml | 42 ++++ recon_sim.py | 34 +++- stack.py | 28 ++- utils.py | 4 + 8 files changed, 757 insertions(+), 3 deletions(-) create mode 100644 cal_instructions.md create mode 100644 cal_response.py create mode 100644 cal_utils.py create mode 100644 input/cal_recon_sim.yml create mode 100644 input/cal_stack.yml diff --git a/cal_instructions.md b/cal_instructions.md new file mode 100644 index 0000000..a120392 --- /dev/null +++ b/cal_instructions.md @@ -0,0 +1,87 @@ +# QE Response Calibration + +`cal_response.py` measures a Fourier-space transfer function +`R(L) = / <|ktrue|^2>` +on paired NFW-lensed Gaussian CMB sims. The same script works for both +`recon_sim.py` (sims) and `stack.py` (data) because all filter +parameters are read from a YAML preset / CLI args. + +## 1. Generate a calibration + +Two presets have been added: + +```bash +# Sim-side calibration (mirrors recon_sim.py constants) +python cal_response.py --config input/cal_recon_sim.yml \ + --output-dir /cal_sim_inp4 \ + --inpaint-radius-arcmin 6.0 + +# Same, no inpainting (for comparison) +python cal_response.py --config input/cal_recon_sim.yml \ + --output-dir /cal_sim_nohole + +# Data-side calibration (mirrors input/defaults.yml + stack.py) +python cal_response.py --config input/cal_stack.yml \ + --output-dir /cal_data_inp4 \ + --inpaint-radius-arcmin 6.0 +``` + +MPI parallel (recommended for ≥500 sims): + +```bash +OMP_NUM_THREADS=2 mpirun -n 16 python cal_response.py \ + --config input/cal_recon_sim.yml \ + --output-dir /cal_sim_inp4 \ + --inpaint-radius-arcmin 6.0 --nsims 2000 +``` + +Any preset key can be overridden on the CLI (flag = key with `_` -> `-`, +e.g. `--hres-lmax 4000`). Common overrides: `--inpaint-radius-arcmin`, +`--mass`, `--nsims`, `--debug` (10 sims). + +### Output + +``` +/cal_*_inp4/ + cal.npz # ell_edges, ell_centers, R_L, R_L_err, auto_power, cross_mean + cal.yaml # full filter parameter set + NFW metadata +``` + +## 2. Apply the calibration + +Pass `--calibration ` to `recon_sim.py` or `stack.py`. The +loader validates that every parameter in `FILTER_KEYS` matches the +run's actual settings; mismatches raise a `ValueError`. + +```bash +# Sim +python recon_sim.py myrun agora halo --cmb-ksz --inpaint 6.0 \ + --calibration /cal_sim_inp4 + +# Mean-field: SAME calibration must be passed +python recon_sim.py myrun agora halo --cmb-ksz --inpaint 6.0 \ + --calibration /cal_sim_inp4 --is-meanfield + +# Data +python stack.py ... --inpaint --calibration /cal_data_inp4 +python stack.py ... --inpaint --calibration /cal_data_inp4 \ + --is-meanfield +``` + +The integration divides `rkmap` by `R_2d` immediately after the symlens +QE reconstruct. Because this is a linear operation, you can subtract the +mean-field at radial-profile time as before — provided both runs used +the same `--calibration` directory. + +## 3. Note + +- **`stack.py` and `recon_sim.py` have different settings.** Do *not* + reuse the same calibration product across them. The validator will + refuse, but it's faster to use the right preset from the start. +- **Foregrounds are intentionally absent in the sims.** `R(L)` depends + on the QE filter weights only, not on the data. Adding foregrounds + would add scatter without changing the mean response. +- **Mean-field with calibration.** Always pass `--calibration` to + *both* the cluster run and the `--is-meanfield` run, or to neither. +- **`R(L)` is mildly mass-dependent.** Calibrate at the median mass of + the sample (default `--mass 3.5e14`). diff --git a/cal_response.py b/cal_response.py new file mode 100644 index 0000000..613a348 --- /dev/null +++ b/cal_response.py @@ -0,0 +1,434 @@ +"""Calibrate the cmbhalolensing QE response on NFW-lensed Gaussian sims. + +Produces a Fourier-space transfer function + R(L) = / <|ktrue_k|^2> +where ``response_k = QE(lensed) - QE(unlensed)`` on paired Gaussian-CMB +sims. The pairing makes the mean field cancel nearly exactly; no foregrounds +or noise are added (R(L) depends only on the QE filter weights, not on +the data itself). + +The same script can produce calibration products for both ``recon_sim.py`` +and ``stack.py`` by changing the filter parameters via ``--config`` or +CLI flags. Two presets are included: ``input/cal_recon_sim.yml`` and +``input/cal_stack.yml``. + +Outputs (under ``--output-dir``): + cal.npz (ell_edges, ell_centers, R_L, R_L_err, auto_power, cross_mean) + cal.yaml (full filter parameter dict + NFW metadata) + +Apply with ``--calibration `` in recon_sim.py / stack.py. + +MPI-parallel via orphics.mpi.distribute. +""" + +import argparse +import os +import time + +import numpy as np +import yaml + +from pixell import enmap, utils, lensing as enlensing +from orphics import cosmology, lensing as olensing, maps, mpi, stats +from symlens.qe import QE + +from inpaint_utils import Inpainter + + +# NFW kappa + +def _gnfw(x): + x = np.asarray(x, dtype=float) + out = np.zeros_like(x) + gt = x > 1 + lt = (x < 1) & (x > 0) + eq = np.abs(x - 1) < 1e-6 + y = x[gt] + out[gt] = (1 / (y**2 - 1)) * ( + 1 - 2 / np.sqrt(y**2 - 1) * np.arctan(np.sqrt((y - 1) / (y + 1))) + ) + y = x[lt] + out[lt] = (1 / (y**2 - 1)) * ( + 1 - 2 / np.sqrt(1 - y**2) * np.arctanh(np.sqrt((1 - y) / (1 + y))) + ) + out[eq] = 1.0 / 3.0 + return out + + +def nfw_kappa_map(shape, wcs, M200c, z, c200c): + """Projected NFW convergence on an enmap grid (colossus cosmology).""" + from colossus.cosmology import cosmology as col_cosmo + col_cosmo.setCosmology("planck18") + cosmo = col_cosmo.getCurrent() + rho_c = cosmo.rho_c(z) + R200c_kpc = (3 * M200c / (4 * np.pi * 200 * rho_c)) ** (1 / 3) + rs_Mpc = R200c_kpc / (c200c * 1000) + chi_l = cosmo.comovingDistance(0, z) + h = cosmo.H0 / 100.0 + chi_s = 14000.0 * h + theta_s = rs_Mpc / chi_l + win = (chi_s - chi_l) / chi_s + fc = np.log(1 + c200c) - c200c / (1 + c200c) + prefactor = 9.571e-20 * chi_l * (1 + z) * win * M200c / (rs_Mpc**2 * fc) + modrmap = enmap.modrmap(shape, wcs) + x = modrmap / theta_s + kappa = enmap.enmap(prefactor * _gnfw(x), wcs) + ny, nx = shape[-2:] + if ny % 2 == 1 and nx % 2 == 1: + cy, cx = ny // 2, nx // 2 + kappa[cy, cx] = kappa[cy - 1, cx] + return kappa + + +# ILC + +def ilc_two_channel(modlmap, k1, k2, p11, p22, p12, b1, b2, ell_min, ell_max): + """Beam-deconvolved ILC of two channels with band restriction.""" + sel = np.logical_and(modlmap >= ell_min, modlmap <= ell_max) + nells = int(sel.sum()) + cov = np.zeros((nells, 2, 2)) + cov[:, 0, 0] = p11[sel] + cov[:, 1, 1] = p22[sel] + cov[:, 0, 1] = p12[sel] + cov[:, 1, 0] = p12[sel] + ms = np.stack([k1[sel], k2[sel]]).swapaxes(0, 1) + rs = np.stack([b1[sel], b2[sel]]).swapaxes(0, 1) + num = np.array([np.linalg.solve(cov[i], ms[i]) for i in range(nells)]) + den = np.array([np.linalg.solve(cov[i], rs[i]) for i in range(nells)]) + tcov = 1.0 / np.einsum("ij,ij->i", rs, den) + out = np.zeros_like(k1) + out[sel] = np.einsum("ij,ij->i", rs, num) * tcov + tret = np.zeros_like(p11) + tret[sel] = tcov + return out, tret + + +# Parameter handling +FILTER_KEYS = ( + "pix_arcmin", "stamp_width_arcmin", + "taper_percent", "pad_percent", + "grad_lmin", "grad_lmax", + "hres_lmin", "hres_lmax", "hres_lxcut", "hres_lycut", + "kappa_lmin", "kappa_lmax", + "grad_beam_fwhm", "grad_noise_uK_arcmin", + "hres_beam_fwhms", "hres_noise_uK_arcmin", + "inpaint_radius_arcmin", + "inpaint_context_arcmin", "inpaint_pcov_res_arcmin", + "qe_estimator", +) + +_DEFAULTS = { + "pix_arcmin": 0.5, + "stamp_width_arcmin": 120.0, + "taper_percent": 20.0, + "pad_percent": 3.0, + "grad_lmin": 200, + "grad_lmax": 2000, + "hres_lmin": 200, + "hres_lmax": 3500, + "hres_lxcut": 2, + "hres_lycut": 2, + "kappa_lmin": 200, + "kappa_lmax": 3000, + "grad_beam_fwhm": 5.0, + "grad_noise_uK_arcmin": 35.0, + "hres_beam_fwhms": [1.5, 2.2], + "hres_noise_uK_arcmin": 15.0, + "inpaint_radius_arcmin": 0.0, + "inpaint_context_arcmin": 60.0, + "inpaint_pcov_res_arcmin": 2.0, + "qe_estimator": "hdv", +} + + +def parse_args(): + ap = argparse.ArgumentParser( + description=__doc__, formatter_class=argparse.RawTextHelpFormatter, + ) + ap.add_argument("--config", default=None, + help="YAML with filter parameters (CLI args override).") + ap.add_argument("--output-dir", required=True) + for k, v in _DEFAULTS.items(): + flag = "--" + k.replace("_", "-") + if isinstance(v, list): + ap.add_argument(flag, type=float, nargs="+", default=None) + elif isinstance(v, bool): + ap.add_argument(flag, action="store_true", default=None) + elif isinstance(v, int) and not isinstance(v, bool): + ap.add_argument(flag, type=int, default=None) + elif isinstance(v, float): + ap.add_argument(flag, type=float, default=None) + else: + ap.add_argument(flag, type=str, default=None) + ap.add_argument("--ell-bin-width", type=int, default=200) + ap.add_argument("--nsims", type=int, default=2000) + ap.add_argument("--mass", type=float, default=3.5e14) + ap.add_argument("--redshift", type=float, default=0.5) + ap.add_argument("--concentration", type=float, default=3.2) + ap.add_argument("--seed", type=int, default=99999) + ap.add_argument("--debug", action="store_true") + return ap.parse_args() + + +def resolve_params(args): + """Merge defaults, --config YAML, and CLI overrides (in that order).""" + p = dict(_DEFAULTS) + if args.config is not None: + with open(args.config) as f: + cfg = yaml.safe_load(f) or {} + for k, v in cfg.items(): + if k in _DEFAULTS: + p[k] = v + for k in _DEFAULTS: + val = getattr(args, k, None) + if val is not None: + p[k] = val + p["ell_bin_width"] = args.ell_bin_width + p["nsims"] = 10 if args.debug else args.nsims + p["mass"] = args.mass + p["redshift"] = args.redshift + p["concentration"] = args.concentration + p["seed"] = args.seed + p["hres_beam_fwhms"] = [float(f) for f in p["hres_beam_fwhms"]] + return p + + +# Sim loop + +def main(): + args = parse_args() + params = resolve_params(args) + + comm, rank, my_tasks = mpi.distribute(params["nsims"]) + size = comm.Get_size() + + if rank == 0: + os.makedirs(args.output_dir, exist_ok=True) + print(f"[cal] output: {args.output_dir}") + print(f"[cal] nsims={params['nsims']}, ranks={size}") + print(f"[cal] stamp={params['stamp_width_arcmin']}', " + f"pix={params['pix_arcmin']}', taper={params['taper_percent']}%") + print(f"[cal] grad: B={params['grad_beam_fwhm']}', " + f"n={params['grad_noise_uK_arcmin']} uK-arcmin, " + f"ell=[{params['grad_lmin']}, {params['grad_lmax']}]") + print(f"[cal] hres: B={params['hres_beam_fwhms']}', " + f"n={params['hres_noise_uK_arcmin']} uK-arcmin, " + f"ell=[{params['hres_lmin']}, {params['hres_lmax']}], " + f"lxcut={params['hres_lxcut']}, lycut={params['hres_lycut']}") + print(f"[cal] inpaint radius: {params['inpaint_radius_arcmin']}'") + + # Geometry + w = (params["stamp_width_arcmin"] / 60.0) / 2 * utils.degree + r = params["pix_arcmin"] * utils.arcmin + shape, wcs = enmap.geometry(pos=[[-w, -w], [w, w]], res=r, proj="car") + modlmap = enmap.modlmap(shape, wcs) + + taper_arr, _ = maps.get_taper(shape, wcs, + taper_percent=params["taper_percent"], + pad_percent=params["pad_percent"], + weight=None) + + # Beams + Bplc = maps.gauss_beam(modlmap, params["grad_beam_fwhm"]) + hres_fwhms = list(params["hres_beam_fwhms"]) + Bhres = [maps.gauss_beam(modlmap, f) for f in hres_fwhms] + n_hres = len(Bhres) + if n_hres not in (1, 2): + raise ValueError(f"hres_beam_fwhms must have length 1 or 2, got {n_hres}") + + # Theory + theory = cosmology.default_theory() + # Match recon_sim.py / stack.py convention: lensed CMB in uC_T_T. + ucltt = theory.lCl("TT", modlmap) + + # Filter masks + xmask = maps.mask_kspace(shape, wcs, + lmin=params["grad_lmin"], lmax=params["grad_lmax"]) + ymask = maps.mask_kspace(shape, wcs, + lmin=params["hres_lmin"], lmax=params["hres_lmax"], + lxcut=params["hres_lxcut"], lycut=params["hres_lycut"]) + kmask = maps.mask_kspace(shape, wcs, + lmin=params["kappa_lmin"], lmax=params["kappa_lmax"]) + + # Total power spectra (theoretical, no foregrounds) + n_plc_pix = (params["grad_noise_uK_arcmin"] * utils.arcmin) ** 2 + tclpp = ucltt + n_plc_pix / Bplc**2 # Planck total, beam-deconvolved + n_hres_pix = (params["hres_noise_uK_arcmin"] * utils.arcmin) ** 2 + + if n_hres == 2: + p11 = ucltt * Bhres[0]**2 + n_hres_pix + p22 = ucltt * Bhres[1]**2 + n_hres_pix + p12 = ucltt * Bhres[0] * Bhres[1] + # Analytic ILC residual var = 1 / (r^T C^-1 r); equals total + # beam-deconvolved post-ILC power. + sel = (modlmap >= params["hres_lmin"]) & (modlmap <= params["hres_lmax"]) + nells = int(sel.sum()) + cov = np.zeros((nells, 2, 2)) + cov[:, 0, 0] = p11[sel]; cov[:, 1, 1] = p22[sel] + cov[:, 0, 1] = p12[sel]; cov[:, 1, 0] = p12[sel] + rs = np.stack([Bhres[0][sel], Bhres[1][sel]]).swapaxes(0, 1) + den = np.array([np.linalg.solve(cov[i], rs[i]) for i in range(nells)]) + tclaa = ucltt.copy() + tclaa[sel] = 1.0 / np.einsum("ij,ij->i", rs, den) + else: + tclaa = ucltt + n_hres_pix / Bhres[0]**2 + + # Inpainter + grad_inpainter = None + if params["inpaint_radius_arcmin"] > 0: + grad_inpainter = Inpainter( + fine_shape=shape, fine_wcs=wcs, + hole_radius_arcmin=params["inpaint_radius_arcmin"], + theory=theory, + beam_fn=lambda x: maps.gauss_beam(x, params["grad_beam_fwhm"]), + noise_uK_arcmin=params["grad_noise_uK_arcmin"], + context_width_arcmin=params["inpaint_context_arcmin"], + pcov_res_arcmin=params["inpaint_pcov_res_arcmin"], + ) + + # NFW kappa + deflection + kappa_true = nfw_kappa_map(shape, wcs, params["mass"], + params["redshift"], params["concentration"]) + grad_phi = olensing.alpha_from_kappa(kappa_true) + ktrue_k = np.asarray( + enmap.fft(enmap.enmap(np.asarray(kappa_true) * taper_arr, wcs), + normalize="phys") * kmask + ) + if rank == 0: + print(f"[cal] NFW: M={params['mass']:.2e}, z={params['redshift']}, " + f"c={params['concentration']}, " + f"peak_kappa={float(np.max(np.asarray(kappa_true))):.4f}") + + # Ell binning + ell_edges = np.arange(params["kappa_lmin"], + params["kappa_lmax"] + params["ell_bin_width"], + params["ell_bin_width"]) + n_ell = len(ell_edges) - 1 + ell_centers = (ell_edges[:-1] + ell_edges[1:]) / 2.0 + ell_binner = stats.bin2D(modlmap, ell_edges) + _, auto_power = ell_binner.bin(np.abs(ktrue_k) ** 2) + + # CMB generator (unlensed) + ps_cmb = np.zeros((1, 1) + tuple(shape)) + ps_cmb[0, 0] = theory.uCl("TT", modlmap) + mgen = maps.MapGen(shape=(1,) + tuple(shape), wcs=wcs, cov=ps_cmb) + + # Sim loop + s = stats.Stats(comm) + t0 = time.time() + n_my = len(my_tasks) + + feed_dict_base = { + "uC_T_T": ucltt, + "tC_A_T_A_T": tclaa, + "tC_P_T_P_T": tclpp, + "tC_A_T_P_T": ucltt, + "tC_P_T_A_T": ucltt, + } + + for j, isim in enumerate(my_tasks): + if rank == 0 and (j % max(1, n_my // 10) == 0): + print(f"[cal] sim {j}/{n_my} ({time.time()-t0:.0f}s)") + seed = params["seed"] + int(isim) + + unlensed = np.asarray(mgen.get_map(seed=seed)[0], dtype=np.float64) + lensed = enlensing.lens_map(unlensed, grad_phi, + order=5, border="cyclic") + + krecon_pair = [] + for cmb_input in (lensed, unlensed): + cmb_in = enmap.enmap(np.asarray(cmb_input, dtype=np.float64), wcs) + + # Gradient leg: beam-convolve, inpaint (untapered), taper, FFT, deconvolve. + g_k = enmap.fft(cmb_in, normalize="phys") * Bplc + g_map = enmap.enmap(enmap.ifft(g_k, normalize="phys").real, wcs) + if grad_inpainter is not None: + g_map = grad_inpainter.inpaint(g_map) + k_grad = (enmap.fft(enmap.enmap(g_map * taper_arr, wcs), + normalize="phys") / Bplc) + + # High-res leg: per-freq beam-convolve, taper, FFT, ILC (or single). + if n_hres == 1: + h_k = enmap.fft(cmb_in, normalize="phys") * Bhres[0] + h_map = enmap.ifft(h_k, normalize="phys").real + k_hres = (enmap.fft(enmap.enmap(h_map * taper_arr, wcs), + normalize="phys") / Bhres[0]) + else: + kfreq = [] + for B in Bhres: + h_k = enmap.fft(cmb_in, normalize="phys") * B + h_map = enmap.ifft(h_k, normalize="phys").real + kfreq.append(np.asarray( + enmap.fft(enmap.enmap(h_map * taper_arr, wcs), + normalize="phys") + )) + k_hres, _ = ilc_two_channel( + modlmap, kfreq[0], kfreq[1], + p11, p22, p12, Bhres[0], Bhres[1], + params["hres_lmin"], params["hres_lmax"], + ) + + feed_dict = dict(feed_dict_base) + feed_dict["X"] = k_grad + feed_dict["Y"] = k_hres + qe_obj = QE(shape, wcs, feed_dict, + estimator=params["qe_estimator"], XY="TT", + xmask=xmask, ymask=ymask, + field_names=["P", "A"], groups=None, kmask=kmask) + rkmap = qe_obj.reconstruct(feed_dict, + xname="X_l1", yname="Y_l2", + physical_units=True) + krecon_pair.append(np.asarray(rkmap)) + + response_k = krecon_pair[0] - krecon_pair[1] + if not np.all(np.isfinite(response_k)): + continue + cross_2d = np.real(response_k * np.conj(ktrue_k)) + _, cross_binned = ell_binner.bin(cross_2d) + s.add_to_stats("cross", cross_binned) + s.add_to_stats("nused", np.array([1.0])) + + s.get_stats() + if rank != 0: + return + + cross_mean = s.stats["cross"]["mean"] + cross_err = s.stats["cross"]["errmean"] + n_used = int(np.round(s.stats["nused"]["mean"] * size)) + + safe = auto_power > 1e-30 + R_L = np.ones(n_ell) + R_L_err = np.zeros(n_ell) + R_L[safe] = cross_mean[safe] / auto_power[safe] + R_L_err[safe] = cross_err[safe] / auto_power[safe] + + # Save + out_yaml = {k: params[k] for k in FILTER_KEYS} + out_yaml["nfw"] = { + "mass": float(params["mass"]), + "redshift": float(params["redshift"]), + "concentration": float(params["concentration"]), + "peak_kappa": float(np.max(np.asarray(kappa_true))), + } + out_yaml["nsims_requested"] = int(params["nsims"]) + out_yaml["nsims_used"] = int(n_used) + out_yaml["ell_bin_width"] = int(params["ell_bin_width"]) + + yaml_path = os.path.join(args.output_dir, "cal.yaml") + with open(yaml_path, "w") as f: + yaml.safe_dump(out_yaml, f, sort_keys=False, default_flow_style=None) + + np.savez( + os.path.join(args.output_dir, "cal.npz"), + ell_edges=ell_edges, ell_centers=ell_centers, + R_L=R_L, R_L_err=R_L_err, + auto_power=auto_power, cross_mean=cross_mean, + ) + print(f"[cal] Saved {yaml_path}") + print(f"[cal] Saved cal.npz ({n_used}/{params['nsims']} sims used)") + print(f"[cal] mean R(L) in [{params['kappa_lmin']},{params['kappa_lmax']}]: " + f"{float(np.mean(R_L[safe])):.4f}") + + +if __name__ == "__main__": + main() diff --git a/cal_utils.py b/cal_utils.py new file mode 100644 index 0000000..eb08e77 --- /dev/null +++ b/cal_utils.py @@ -0,0 +1,98 @@ +"""Load, validate, and apply QE response calibration products. + +Used by ``recon_sim.py`` and ``stack.py`` when ``--calibration `` is +passed. The directory must contain ``cal.npz`` and ``cal.yaml`` written +by ``cal_response.py``. +""" + +import os + +import numpy as np +import yaml + + +FILTER_KEYS = ( + "pix_arcmin", "stamp_width_arcmin", + "taper_percent", "pad_percent", + "grad_lmin", "grad_lmax", + "hres_lmin", "hres_lmax", "hres_lxcut", "hres_lycut", + "kappa_lmin", "kappa_lmax", + "grad_beam_fwhm", "grad_noise_uK_arcmin", + "hres_beam_fwhms", "hres_noise_uK_arcmin", + "inpaint_radius_arcmin", + "inpaint_context_arcmin", "inpaint_pcov_res_arcmin", + "qe_estimator", +) + + +def load_calibration(cal_dir): + """Return {'ell_centers', 'R_L', 'R_L_err', 'params', 'path'}.""" + npz_path = os.path.join(cal_dir, "cal.npz") + yaml_path = os.path.join(cal_dir, "cal.yaml") + if not (os.path.isfile(npz_path) and os.path.isfile(yaml_path)): + raise FileNotFoundError( + f"Calibration directory missing cal.npz or cal.yaml: {cal_dir}" + ) + data = np.load(npz_path) + with open(yaml_path) as f: + params = yaml.safe_load(f) + return { + "ell_centers": np.asarray(data["ell_centers"]), + "R_L": np.asarray(data["R_L"]), + "R_L_err": np.asarray(data["R_L_err"]), + "params": params, + "path": os.path.abspath(cal_dir), + } + + +def _equal(a, b): + if a is None or b is None: + return a is b + if isinstance(a, (list, tuple)) or isinstance(b, (list, tuple)): + try: + a = [float(v) for v in a] + b = [float(v) for v in b] + except (TypeError, ValueError): + return list(a) == list(b) + if len(a) != len(b): + return False + return all(np.isclose(x, y, rtol=1e-6, atol=1e-9) for x, y in zip(a, b)) + if isinstance(a, float) or isinstance(b, float): + return np.isclose(float(a), float(b), rtol=1e-6, atol=1e-9) + return a == b + + +def check_compatibility(cal_params, **run): + """Raise ValueError if any filter parameter in run mismatches calibration.""" + mismatches = [] + for key in FILTER_KEYS: + if key not in run: + continue + a = cal_params.get(key) + b = run[key] + if not _equal(a, b): + mismatches.append(f" {key}: calibration={a!r}, run={b!r}") + if mismatches: + raise ValueError( + "Calibration is incompatible with this run:\n" + "\n".join(mismatches) + ) + + +def build_R_2d(R_L, ell_centers, modlmap, fill_value=1.0): + """Linear interpolation of R(L) onto the 2D Fourier grid. + + Outside [ell_centers[0], ell_centers[-1]] -> ``fill_value`` (default 1.0, + a no-op when dividing). Bands killed by kmask are unaffected. + """ + return np.interp( + np.asarray(modlmap), + np.asarray(ell_centers), + np.asarray(R_L), + left=fill_value, + right=fill_value, + ) + + +def apply_R_2d(rkmap, R_2d, eps=1e-3): + """Safe division: rkmap / R_2d, zero where |R_2d| < eps.""" + return np.where(np.abs(R_2d) > eps, rkmap / R_2d, 0.0) diff --git a/input/cal_recon_sim.yml b/input/cal_recon_sim.yml new file mode 100644 index 0000000..73a350b --- /dev/null +++ b/input/cal_recon_sim.yml @@ -0,0 +1,33 @@ +# Calibration preset that mirrors recon_sim.py's hardcoded settings. +# Pass via: python cal_response.py --config input/cal_recon_sim.yml ... +# Override the inpaint radius on the CLI when calibrating different holes. + +pix_arcmin: 0.5 +stamp_width_arcmin: 120.0 + +taper_percent: 20.0 +pad_percent: 3.0 + +grad_lmin: 200 +grad_lmax: 2000 + +hres_lmin: 200 +hres_lmax: 3500 +hres_lxcut: 2 +hres_lycut: 2 + +kappa_lmin: 200 +kappa_lmax: 3000 + +grad_beam_fwhm: 5.0 +grad_noise_uK_arcmin: 35.0 + +# Two-channel ILC; matches ilc_beam_fwhm150 / ilc_beam_fwhm090 in recon_sim.py. +hres_beam_fwhms: [1.5, 2.2] +hres_noise_uK_arcmin: 15.0 + +inpaint_radius_arcmin: 0.0 # override on CLI (e.g. --inpaint-radius-arcmin 6.0) +inpaint_context_arcmin: 60.0 +inpaint_pcov_res_arcmin: 2.0 + +qe_estimator: hdv diff --git a/input/cal_stack.yml b/input/cal_stack.yml new file mode 100644 index 0000000..dc25f4a --- /dev/null +++ b/input/cal_stack.yml @@ -0,0 +1,42 @@ +# Calibration preset that mirrors input/defaults.yml + stack.py setup. +# Pass via: python cal_response.py --config input/cal_stack.yml ... +# Override the inpaint radius on the CLI when calibrating different holes. + +pix_arcmin: 0.5 +stamp_width_arcmin: 128.0 + +taper_percent: 12.0 +pad_percent: 3.0 + +# From defaults.yml: gradient_lmin / gradient_lmax +grad_lmin: 200 +grad_lmax: 2000 + +# From defaults.yml: conservative_highres_* (set to aggressive_* if running +# stack.py without --conservative; calibrate one per filter). +hres_lmin: 500 +hres_lmax: 6000 +hres_lxcut: 50 +hres_lycut: 90 + +# From defaults.yml: kappa_Lmin / kappa_Lmax +kappa_lmin: 200 +kappa_lmax: 5000 + +# planck_smica_beam_fwhm: 5.0; gradient_fiducial_rms 40 in defaults.yml -- +# recon_sim.py uses 35; pick what your stack.py actually runs with. +grad_beam_fwhm: 5.0 +grad_noise_uK_arcmin: 35.0 + +# stack.py loads real per-frequency beams from disk (cutils.load_beam), +# but for the response calibration a Gaussian approximation is adequate +# (the QE filter weights are insensitive to fine beam shape over the +# ILC ell band). Change here if you have a better effective fwhm. +hres_beam_fwhms: [1.5, 2.2] +hres_noise_uK_arcmin: 20.0 # defaults.yml: highres_fiducial_rms + +inpaint_radius_arcmin: 0.0 # override on CLI (e.g. --inpaint-radius-arcmin 6.0) +inpaint_context_arcmin: 60.0 +inpaint_pcov_res_arcmin: 2.0 + +qe_estimator: hdv diff --git a/recon_sim.py b/recon_sim.py index f7dca0d..b9d5eb8 100644 --- a/recon_sim.py +++ b/recon_sim.py @@ -86,7 +86,11 @@ parser.add_argument( "--inpaint", type=float, default=None, help="Inpainting for gradient leg. Enter the inpainting radius in arcmin." ) -args = parser.parse_args() +parser.add_argument( + "--calibration", type=str, default=None, + help="Directory written by cal_response.py. Divides krecon by R(L) before stacking." +) +args = parser.parse_args() if args.which_sim == "websky": output_path = paths.websky_output_path elif args.which_sim == "sehgal": output_path = paths.sehgal_output_path @@ -1028,10 +1032,32 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): l_edges = np.arange(minell / 2, 8001, minell) lbinner = stats.bin2D(modlmap, l_edges) # PS correction factor - w2 = np.mean(taper ** 2) + w2 = np.mean(taper ** 2) # fsky for bandpower variance fsky = enmap.area(shape, wcs) * w2 / 4.0 / np.pi + R_2d = None + if args.calibration is not None: + from cal_utils import load_calibration, check_compatibility, build_R_2d + _cal = load_calibration(args.calibration) + check_compatibility( + _cal["params"], + pix_arcmin=px, + stamp_width_arcmin=width_deg * 60, + taper_percent=tap_per, pad_percent=pad_per, + grad_lmin=xlmin, grad_lmax=xlmax, + hres_lmin=ylmin, hres_lmax=ylmax, + hres_lxcut=lycut, hres_lycut=lycut, + kappa_lmin=klmin, kappa_lmax=klmax, + grad_beam_fwhm=fwhm_plc, grad_noise_uK_arcmin=nlevel_plc, + hres_beam_fwhms=[ilc_beam_fwhm150, ilc_beam_fwhm090], + hres_noise_uK_arcmin=nlevel, + inpaint_radius_arcmin=(args.inpaint or 0.0), + ) + R_2d = build_R_2d(_cal["R_L"], _cal["ell_centers"], modlmap) + if rank == 0: + print(f" ::: applying response calibration from {_cal['path']}") + # same filter as the post-reconstuction for true kappa (also tapered!) k_stamp = maps.filter_map(kstamp*taper, kmask) s.add_to_stack("kstamp", k_stamp) @@ -1293,6 +1319,10 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): # Fourier space lens reconstruction rkmap = cqe.reconstruct(feed_dict, xname="X_l1", yname="Y_l2", physical_units=True) + if R_2d is not None: + from cal_utils import apply_R_2d + rkmap = apply_R_2d(rkmap, R_2d) + assert np.all(np.isfinite(rkmap)) diff --git a/stack.py b/stack.py index f8c9e4f..07d6990 100644 --- a/stack.py +++ b/stack.py @@ -879,6 +879,28 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): # map of distances from center modrmap = enmap.modrmap(shape, wcs) + R_2d = None + if args.calibration is not None: + from cal_utils import load_calibration, check_compatibility, build_R_2d + _cal = load_calibration(args.calibration) + check_compatibility( + _cal["params"], + pix_arcmin=args.pwidth, + stamp_width_arcmin=args.swidth, + taper_percent=args.tap_per, pad_percent=args.pad_per, + grad_lmin=xlmin, grad_lmax=xlmax, + hres_lmin=ylmin, hres_lmax=ylmax, + hres_lxcut=args.hres_lxcut, hres_lycut=args.hres_lycut, + kappa_lmin=klmin, kappa_lmax=klmax, + grad_beam_fwhm=plc_beam_fwhm, + grad_noise_uK_arcmin=defaults.gradient_fiducial_rms, + hres_noise_uK_arcmin=defaults.highres_fiducial_rms, + inpaint_radius_arcmin=(4.0 if args.inpaint else 0.0), + ) + R_2d = build_R_2d(_cal["R_L"], _cal["ell_centers"], modlmap) + if rank == 0: + print(f" ::: applying response calibration from {_cal['path']}") + # Fourier map -> PS pow = lambda x, y: (x * y.conj()).real @@ -1118,7 +1140,11 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): ) # Fourier space lens reconstruction krecon = cqe.reconstruct(feed_dict, xname="X_l1", yname="Y_l2", physical_units=True) - # not the pixel unit + # not the pixel unit + + if R_2d is not None: + from cal_utils import apply_R_2d + krecon = apply_R_2d(krecon, R_2d) # transform to real space for unweighted stack kappa = enmap.ifft(krecon, normalize="phys").real diff --git a/utils.py b/utils.py index fddb13c..abe43ad 100644 --- a/utils.py +++ b/utils.py @@ -154,6 +154,10 @@ def initialize_pipeline_config(): parser.add_argument( "--inpaint", action="store_true", help="Inpaint gradient." ) + parser.add_argument( + "--calibration", type=str, default=None, + help="Directory written by cal_response.py. Divides krecon by R(L) before stacking." + ) parser.add_argument( "--no-sz-sub", action="store_true", From bc7e0eab299a5e349c6c802837f66a51e54daf3a Mon Sep 17 00:00:00 2001 From: Mathew Madhavacheril Date: Tue, 19 May 2026 23:39:39 -0400 Subject: [PATCH 07/14] cal fixes --- cal_response.py | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/cal_response.py b/cal_response.py index 613a348..ed616ed 100644 --- a/cal_response.py +++ b/cal_response.py @@ -58,7 +58,7 @@ def _gnfw(x): def nfw_kappa_map(shape, wcs, M200c, z, c200c): """Projected NFW convergence on an enmap grid (colossus cosmology).""" from colossus.cosmology import cosmology as col_cosmo - col_cosmo.setCosmology("planck18") + col_cosmo.setCosmology("planck18",persistence="") cosmo = col_cosmo.getCurrent() rho_c = cosmo.rho_c(z) R200c_kpc = (3 * M200c / (4 * np.pi * 200 * rho_c)) ** (1 / 3) @@ -290,10 +290,8 @@ def main(): kappa_true = nfw_kappa_map(shape, wcs, params["mass"], params["redshift"], params["concentration"]) grad_phi = olensing.alpha_from_kappa(kappa_true) - ktrue_k = np.asarray( - enmap.fft(enmap.enmap(np.asarray(kappa_true) * taper_arr, wcs), - normalize="phys") * kmask - ) + ktrue_k = enmap.fft(enmap.enmap(kappa_true * taper_arr, wcs),normalize="phys") * kmask + if rank == 0: print(f"[cal] NFW: M={params['mass']:.2e}, z={params['redshift']}, " f"c={params['concentration']}, " @@ -331,13 +329,12 @@ def main(): print(f"[cal] sim {j}/{n_my} ({time.time()-t0:.0f}s)") seed = params["seed"] + int(isim) - unlensed = np.asarray(mgen.get_map(seed=seed)[0], dtype=np.float64) + unlensed = mgen.get_map(seed=seed)[0] lensed = enlensing.lens_map(unlensed, grad_phi, order=5, border="cyclic") krecon_pair = [] - for cmb_input in (lensed, unlensed): - cmb_in = enmap.enmap(np.asarray(cmb_input, dtype=np.float64), wcs) + for cmb_in in (lensed, unlensed): # Gradient leg: beam-convolve, inpaint (untapered), taper, FFT, deconvolve. g_k = enmap.fft(cmb_in, normalize="phys") * Bplc @@ -378,7 +375,7 @@ def main(): rkmap = qe_obj.reconstruct(feed_dict, xname="X_l1", yname="Y_l2", physical_units=True) - krecon_pair.append(np.asarray(rkmap)) + krecon_pair.append(rkmap) response_k = krecon_pair[0] - krecon_pair[1] if not np.all(np.isfinite(response_k)): @@ -394,7 +391,7 @@ def main(): cross_mean = s.stats["cross"]["mean"] cross_err = s.stats["cross"]["errmean"] - n_used = int(np.round(s.stats["nused"]["mean"] * size)) + n_used = int(np.round(float(s.stats["nused"]["mean"][0]) * size)) safe = auto_power > 1e-30 R_L = np.ones(n_ell) From 36e9f24fd871448dde41afa0f87cdb061c084861 Mon Sep 17 00:00:00 2001 From: Mathew Madhavacheril Date: Tue, 19 May 2026 23:42:25 -0400 Subject: [PATCH 08/14] calibration response plots --- cal_response.py | 133 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 133 insertions(+) diff --git a/cal_response.py b/cal_response.py index ed616ed..2020852 100644 --- a/cal_response.py +++ b/cal_response.py @@ -385,7 +385,14 @@ def main(): s.add_to_stats("cross", cross_binned) s.add_to_stats("nused", np.array([1.0])) + # Real-space response stamp (for diagnostic plots only) + response_real = np.asarray( + enmap.ifft(enmap.enmap(response_k, wcs), normalize="phys").real + ) + s.add_to_stack("response_real", response_real) + s.get_stats() + s.get_stacks() if rank != 0: return @@ -426,6 +433,132 @@ def main(): print(f"[cal] mean R(L) in [{params['kappa_lmin']},{params['kappa_lmax']}]: " f"{float(np.mean(R_L[safe])):.4f}") + # Diagnostic plots + try: + response_real_stack = np.asarray(s.stacks["response_real"]) + kappa_filt_real = np.asarray( + enmap.ifft(enmap.enmap(ktrue_k, wcs), normalize="phys").real + ) + modrmap_arcmin = enmap.modrmap(shape, wcs) * (180.0 * 60.0 / np.pi) + _make_plots( + args.output_dir, params, ell_centers, R_L, R_L_err, + response_real_stack, kappa_filt_real, modrmap_arcmin, + n_used, + ) + print(f"[cal] Saved plots to {args.output_dir}") + except Exception as e: + print(f"[cal] WARNING: plot generation failed: {e}") + + +def _make_plots(out_dir, params, ell_centers, R_L, R_L_err, + response_real, kappa_filt, modrmap_arcmin, n_used): + """Three diagnostic plots written under out_dir. + + cal_R_L.png : R(L) +/- err vs L, with unity reference + cal_stamps.png : zoomed 2D stamps of filtered true kappa vs mean + recovered response + cal_profile.png : 1D radial profiles of the same + """ + import matplotlib + matplotlib.use("Agg") + import matplotlib.pyplot as plt + + title_suffix = ( + f"hole={params['inpaint_radius_arcmin']}', " + f"M={params['mass']:.2e}, z={params['redshift']}, " + f"nsims={n_used}" + ) + + # ---- R(L) ---- + fig, ax = plt.subplots(figsize=(7.5, 4.5)) + ax.errorbar(ell_centers, R_L, yerr=R_L_err, fmt="o-", color="#1f77b4", + capsize=2, lw=1.5, ms=5, label=r"$R(L)$") + ax.axhline(1.0, ls="--", color="k", alpha=0.5, lw=1, label="unity") + safe = R_L_err < 1.0 # avoid garbage bins outside the kmask band + if safe.any(): + mean_R = float(np.mean(R_L[safe])) + ax.axhline(mean_R, ls=":", color="#d62728", lw=1.5, + label=fr"$\langle R \rangle = {mean_R:.3f}$") + ax.set_xlabel(r"multipole $L$") + ax.set_ylabel(r"$R(L)$") + ax.set_title("QE response transfer function\n" + title_suffix, + fontsize=10) + ax.legend(loc="best", framealpha=0.9) + ax.grid(alpha=0.3) + fig.tight_layout() + fig.savefig(os.path.join(out_dir, "cal_R_L.png"), dpi=130) + plt.close(fig) + + # ---- 2D stamps (zoomed) ---- + n = response_real.shape[-1] + half_arcmin = 20.0 + half_pix = int(half_arcmin / params["pix_arcmin"]) + cy = n // 2 + sl = (slice(cy - half_pix, cy + half_pix), + slice(cy - half_pix, cy + half_pix)) + extent = [-half_arcmin, half_arcmin, -half_arcmin, half_arcmin] + + true_crop = kappa_filt[sl] + resp_crop = response_real[sl] + vmax = max(np.nanmax(true_crop), np.nanmax(resp_crop)) + vmin = min(np.nanmin(true_crop), np.nanmin(resp_crop)) + abs_max = max(abs(vmin), abs(vmax)) + norm_kwargs = dict(vmin=-abs_max, vmax=abs_max, cmap="RdBu_r") + + fig, axs = plt.subplots(1, 3, figsize=(13.5, 4.4), + gridspec_kw={"width_ratios": [1, 1, 1]}) + for ax_, data, ttl in zip( + axs, + [true_crop, resp_crop, resp_crop - true_crop], + ["filtered true kappa", "mean recovered response", + "residual (recov - true)"], + ): + im = ax_.imshow(data, origin="lower", extent=extent, **norm_kwargs) + ax_.set_title(ttl, fontsize=10) + ax_.set_xlabel("arcmin") + plt.colorbar(im, ax=ax_, fraction=0.046) + axs[0].set_ylabel("arcmin") + fig.suptitle(title_suffix, fontsize=10) + fig.tight_layout(rect=(0, 0, 1, 0.95)) + fig.savefig(os.path.join(out_dir, "cal_stamps.png"), dpi=130) + plt.close(fig) + + # ---- 1D radial profile ---- + from orphics.stats import bin2D + r_edges = np.arange(0.0, 15.01, 1.0) + r_centers = (r_edges[:-1] + r_edges[1:]) / 2.0 + rbinner = bin2D(modrmap_arcmin, r_edges) + _, prof_true = rbinner.bin(kappa_filt) + _, prof_resp = rbinner.bin(response_real) + + fig, axs = plt.subplots(2, 1, figsize=(7.5, 6.5), + gridspec_kw={"height_ratios": [3, 1.4]}, + sharex=True) + axs[0].plot(r_centers, prof_true, "o-", color="k", lw=1.5, ms=5, + label="filtered true kappa") + axs[0].plot(r_centers, prof_resp, "s-", color="#1f77b4", lw=1.5, ms=5, + label="mean recovered response") + axs[0].axhline(0, ls="-", color="k", alpha=0.2, lw=1) + axs[0].set_ylabel(r"$\kappa(\theta)$ profile") + axs[0].set_title("Real-space response vs true kappa\n" + title_suffix, + fontsize=10) + axs[0].legend(loc="best", framealpha=0.9) + axs[0].grid(alpha=0.3) + + safe_r = np.abs(prof_true) > 1e-10 + ratio = np.full_like(prof_true, np.nan) + ratio[safe_r] = prof_resp[safe_r] / prof_true[safe_r] + axs[1].plot(r_centers, ratio, "o-", color="#d62728", lw=1.5, ms=5) + axs[1].axhline(1.0, ls="--", color="k", alpha=0.5, lw=1) + axs[1].set_xlabel(r"$\theta$ (arcmin)") + axs[1].set_ylabel("recov / true") + axs[1].set_ylim(0.0, 1.5) + axs[1].grid(alpha=0.3) + + fig.tight_layout() + fig.savefig(os.path.join(out_dir, "cal_profile.png"), dpi=130) + plt.close(fig) + if __name__ == "__main__": main() From 4423c8e713f14e4b61d95c57bd9f41cd7a9f8035 Mon Sep 17 00:00:00 2001 From: Mathew Madhavacheril Date: Sat, 23 May 2026 14:25:30 -0400 Subject: [PATCH 09/14] fix wcs --- cal_utils.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/cal_utils.py b/cal_utils.py index eb08e77..69667a6 100644 --- a/cal_utils.py +++ b/cal_utils.py @@ -94,5 +94,11 @@ def build_R_2d(R_L, ell_centers, modlmap, fill_value=1.0): def apply_R_2d(rkmap, R_2d, eps=1e-3): - """Safe division: rkmap / R_2d, zero where |R_2d| < eps.""" - return np.where(np.abs(R_2d) > eps, rkmap / R_2d, 0.0) + """Safe division: rkmap / R_2d, zero where |R_2d| < eps. + + Preserves the enmap subclass (and wcs) of ``rkmap`` by avoiding + np.where, which strips the subclass. Bands where |R_2d| < eps get + R_safe = inf, so rkmap / R_safe = 0 there. + """ + R_safe = np.where(np.abs(R_2d) > eps, R_2d, np.inf) + return rkmap / R_safe From f501c0b6a647a8aaf3c75002d6090105f3dba4c8 Mon Sep 17 00:00:00 2001 From: Mathew Madhavacheril Date: Wed, 10 Jun 2026 00:11:12 -0400 Subject: [PATCH 10/14] move stamp exclusion to before statistics --- recon_sim.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/recon_sim.py b/recon_sim.py index b9d5eb8..6addaa0 100644 --- a/recon_sim.py +++ b/recon_sim.py @@ -1326,6 +1326,15 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): assert np.all(np.isfinite(rkmap)) + # transform to real space + kappa = enmap.ifft(rkmap, normalize="phys").real + + # filter weird stamps (temporary2) + if np.any(np.abs(kappa) > 15): + print(f"{task} has large kappa") + continue + + # for cross-correlation check ---------------------------------------------- ikmap = enmap.fft(k_stamp, normalize="phys") # k_stamp is tapered and filtered true kappa @@ -1376,13 +1385,6 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): - # transform to real space - kappa = enmap.ifft(rkmap, normalize="phys").real - - # filter weird stamps (temporary2) - if np.any(np.abs(kappa) > 15): - print(f"{task} has large kappa") - continue # stack reconstructed kappa From 3960e0472869681f0096a1c3d832c3fb7513c329 Mon Sep 17 00:00:00 2001 From: Mathew Madhavacheril Date: Wed, 10 Jun 2026 00:14:24 -0400 Subject: [PATCH 11/14] move statistics --- recon_sim.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/recon_sim.py b/recon_sim.py index 6addaa0..cf56616 100644 --- a/recon_sim.py +++ b/recon_sim.py @@ -1060,9 +1060,7 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): # same filter as the post-reconstuction for true kappa (also tapered!) k_stamp = maps.filter_map(kstamp*taper, kmask) - s.add_to_stack("kstamp", k_stamp) binned_true = bin(k_stamp, modrmap * (180 * 60 / np.pi), bin_edges) - s.add_to_stats("tk1d", binned_true) if args.cmb: hres = h_cmb @@ -1390,7 +1388,10 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): # stack reconstructed kappa s.add_to_stack("lstamp", kappa) binned_kappa = bin(kappa, modrmap * (180 * 60 / np.pi), bin_edges) - s.add_to_stats("k1d", binned_kappa) + s.add_to_stats("k1d", binned_kappa) + s.add_to_stack("kstamp", k_stamp) + s.add_to_stats("tk1d", binned_true) + # stack stamps pre-reconstruction as well (same filter as hres leg) if args.cmb: From e0e4a43868f26df5b10b333fb7b4cee0a8d753c7 Mon Sep 17 00:00:00 2001 From: Mathew Madhavacheril Date: Wed, 10 Jun 2026 00:18:31 -0400 Subject: [PATCH 12/14] move exclusion --- stack.py | 48 ++++++++++++++++++++++++------------------------ 1 file changed, 24 insertions(+), 24 deletions(-) diff --git a/stack.py b/stack.py index 07d6990..fc5cf08 100644 --- a/stack.py +++ b/stack.py @@ -1149,6 +1149,30 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): # transform to real space for unweighted stack kappa = enmap.ifft(krecon, normalize="phys").real + """ + !! REJECT WEIRD KAPPA + """ + + ##### temporary 3: to get rid of stamps with tSZ cluster in random locations + if np.any(np.abs(kappa) > 15): + print(f"{task} has large kappa") + if args.debug_anomalies: + io.plot_img( + kappa, + f"{paths.debugdir}kappa_err_stamp_large_kappa_{task}.png", + arc_width=args.swidth, + xlabel="$\\theta_x$ (arcmin)", + ylabel="$\\theta_y$ (arcmin)", + ) + if args.debug_anomalies: + io.plot_img( + astamp_150, + f"{paths.debugdir}act_150_err_stamp_large_kappa_{task}.png", + arc_width=args.swidth, + xlabel="$\\theta_x$ (arcmin)", + ylabel="$\\theta_y$ (arcmin)", + ) + continue @@ -1192,30 +1216,6 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): #--------------------------------------------------------------------------- - """ - !! REJECT WEIRD KAPPA - """ - - ##### temporary 3: to get rid of stamps with tSZ cluster in random locations - if np.any(np.abs(kappa) > 15): - print(f"{task} has large kappa") - if args.debug_anomalies: - io.plot_img( - kappa, - f"{paths.debugdir}kappa_err_stamp_large_kappa_{task}.png", - arc_width=args.swidth, - xlabel="$\\theta_x$ (arcmin)", - ylabel="$\\theta_y$ (arcmin)", - ) - if args.debug_anomalies: - io.plot_img( - astamp_150, - f"{paths.debugdir}act_150_err_stamp_large_kappa_{task}.png", - arc_width=args.swidth, - xlabel="$\\theta_x$ (arcmin)", - ylabel="$\\theta_y$ (arcmin)", - ) - continue # Save power spectra to look at later if args.debug_powers: From 2b84702aa28f1f89b7b61c35bcc6218e18dfcace Mon Sep 17 00:00:00 2001 From: Eunseong Lee <36232273+eunseongleee@users.noreply.github.com> Date: Thu, 11 Jun 2026 22:08:53 -0400 Subject: [PATCH 13/14] check point current version of the scripts used for model subtracted, inpainted, calibrated run for CMB+tSZ+kSZ sample --- here_mcut_agora.py | 69 +++++++++ here_merge_agora.py | 94 ++++++++++++ here_prepare_agora.py | 216 ++++++++++++++++++++++++++ input/sim_data.yml | 73 +++++---- make_sim_suite_act.py | 348 ++++++++++++++++++++++++++++++++++++++++++ recon_sim.py | 10 +- stack.py | 31 ++-- tsz_subtract.py | 10 +- 8 files changed, 795 insertions(+), 56 deletions(-) create mode 100644 here_mcut_agora.py create mode 100644 here_merge_agora.py create mode 100644 here_prepare_agora.py create mode 100644 make_sim_suite_act.py diff --git a/here_mcut_agora.py b/here_mcut_agora.py new file mode 100644 index 0000000..8deb3b0 --- /dev/null +++ b/here_mcut_agora.py @@ -0,0 +1,69 @@ +import numpy as np +import pandas as pd +from astropy.table import Table + + +halo_path = "/data5/sims/agora_sims/agora_halo_forNEMO_partial.parquet.gzip" +lens_path = "/data5/sims/agora_sims/agora_lensedhalo_forNEMO.parquet.gzip" +output_path = "data/agora_halo_1e13cut.fits" +# output_path = "data/agora_halo_5e13cut.fits" + + +H0 = 67.77 +min_mass = 1e13 # Msun + +print("Reading Agora halo catalog...") + +halo_cat = pd.read_parquet(halo_path) +lens_cat = pd.read_parquet(lens_path) + + +zs = halo_cat["totz"].to_numpy() +RAs = halo_cat["totra"].to_numpy() +decs = halo_cat["totdec"].to_numpy() + +M200c = halo_cat["totm200"].to_numpy() # Msun/h +M500c = halo_cat["totm500"].to_numpy() # Msun/h + +RAs_lens = lens_cat["lenra"].to_numpy() +decs_lens = lens_cat["lendec"].to_numpy() + + +# Unit conversion (Msun/h -> Msun) +M200c = M200c / (0.01 * H0) +M500c = M500c / (0.01 * H0) + +print("Initial sample size:", len(zs)) + + +mask = M200c > min_mass + +RAs = RAs[mask] +decs = decs[mask] +zs = zs[mask] +M200c = M200c[mask] +M500c = M500c[mask] + +RAs_lens = RAs_lens[mask] +decs_lens = decs_lens[mask] + +print("Sample size after mass cut:", len(zs), min_mass/1e14) + + +out_tab = Table() + +out_tab["RADeg_unl"] = RAs +out_tab["decDeg_unl"] = decs +out_tab["RADeg"] = RAs_lens +out_tab["decDeg"] = decs_lens +out_tab["redshift"] = zs +out_tab["true_M200c"] = M200c / 1e14 +out_tab["true_M500c"] = M500c / 1e14 + + +print(f"Writing catalogue to {output_path}") +out_tab.write(output_path, overwrite=True) + +print("Done.") + + diff --git a/here_merge_agora.py b/here_merge_agora.py new file mode 100644 index 0000000..137d345 --- /dev/null +++ b/here_merge_agora.py @@ -0,0 +1,94 @@ +import numpy as np +import pandas as pd +import time as t + +start = t.time() + +my_path = "/data5/sims/agora_sims/" + +halo_path = "/data5/sims/agora_sims/full/halocat/" +lens_path = "/home3/eunseong/lens/" + +""" +Slice number ranges from 4-200 and each integer increases by 25 Mpc/h +In other words: +X=0 -> 0-25 Mpc/h +X=1 -> 25-100 Mpc/h + +With lensed coordinates: +Slice number ranges from 4-177 +columns are lensed RA / lensed DEC / magnification / real / image part of rotation +the length should match with the unlensed halo cat per slice +so you can extract masses and other info by merging the two +""" + +total_num = 174 + +halo_dfs = [] +lens_dfs = [] + +for slice_num in range(4, 4+total_num): + + print(f"Processing slice {slice_num}", flush=True) + + npz = np.load(halo_path + f"agora_halolc_rot_{slice_num}_v050223.npz") + + halo_df = pd.DataFrame.from_dict( + {item: npz[item] for item in npz.files} + ) + + halo_dfs.append( + pd.DataFrame( + halo_df, + columns=[ + "totz", + "totra", + "totdec", + "totm200", # x200 critial density [Msun/h] + "totm500", # x500 critial density [Msun/h] + ], + dtype="float32", + ) + ) + + npz.close() + + + npy = np.load(lens_path + f"agora_halos_lenra_lendec_mag_rotreal_rotimag_deflectnside16384_{slice_num}.npy") + + lens_dfs.append( + pd.DataFrame( + npy[:, 0:2], + columns=["lenra", "lendec"], + dtype="float32", + ) + ) + + del npy + + + + +print("Concatenating halo catalog...", flush=True) + +halo_cat = pd.concat(halo_dfs, ignore_index=True) +print(halo_cat, flush=True) + +halo_out = my_path + "agora_halo_forNEMO_partial.parquet.gzip" +halo_cat.to_parquet(halo_out, compression="gzip") + +print(f"Saved: {halo_out}", flush=True) + + +print("Concatenating lensed halo catalog...", flush=True) + +lens_cat = pd.concat(lens_dfs, ignore_index=True) +print(lens_cat, flush=True) + +lens_out = my_path + "agora_lensedhalo_forNEMO.parquet.gzip" +lens_cat.to_parquet(lens_out, compression="gzip") + +print(f"Saved: {lens_out}", flush=True) + +elapsed = t.time() - start +print(f"\n::: reading and merging catalogues took {elapsed:.1f} seconds", flush=True) \ No newline at end of file diff --git a/here_prepare_agora.py b/here_prepare_agora.py new file mode 100644 index 0000000..bc56f30 --- /dev/null +++ b/here_prepare_agora.py @@ -0,0 +1,216 @@ +import numpy as np +from pixell import enmap, utils, bunch +import utils as cutils +from orphics import io, catalogs +from astropy.io import fits + + +paths = cutils.paths +sim_paths = bunch.Bunch(io.config_from_yaml("input/sim_data.yml")) + +sim_version = "nemo-sim-kit61" +sim_type = "agora" + +output_path = f"{sim_paths.cat_path}/{sim_version}/" +io.mkdir(output_path) +print(" ::: saving to", output_path) + + + +def load_cat(file_name): + print(" ::: input tSZ cat:", file_name) + with fits.open(file_name) as hdu: + data = hdu[1].data + + catalog = { + "zs": data["redshift"], + "ras": data["RADeg"], + "decs": data["decDeg"], + "true_M200c": data["true_M200c"], + "true_M500c": data["true_M500c"], + "M200c": data["M200c"], + "snrs": data["SNR"], + } + + io.save_cols( + output_path + f"{sim_type}_tsz_true.txt", + (catalog["ras"], catalog["decs"], catalog["zs"], catalog["true_M200c"]), + ) + + print(f"\nN={len(catalog['zs'])} [min SNR={catalog['snrs'].min():.1f}]") + print(f"mean true M200c={catalog['true_M200c'].mean():.2f}") + print(f"mean true M500c={catalog['true_M500c'].mean():.2f}") + print(f"mean true z={catalog['zs'].mean():.2f}\n") + + return catalog + + + +def generate_random_cat(n_per_obj, ras, choice="halo"): + + seed = 100 + np.random.seed(seed) + + N = n_per_obj * len(ras) + print(f" ::: total number of random objects = {n_per_obj} x {len(ras)} = {N} [{choice}]") + + if choice == "halo": + + dec_min=-np.pi/2 + dec_max=np.pi/2 + ra_min=0 + ra_max=2.*np.pi + + coords = np.zeros((2,N)) + dmin = np.cos(np.pi/2 - dec_min) + dmax = np.cos(np.pi/2 - dec_max) + coords[0,:] = (np.pi/2 - np.arccos(np.random.uniform(dmin, dmax, N))) / utils.degree + coords[1,:] = np.random.uniform(ra_min, ra_max, N) / utils.degree + + elif choice == "tsz": + + mask_path = f"{paths.data}mask/websky_act_mask.fits" + mask = enmap.read_map(mask_path) + shape, wcs = mask.shape, mask.wcs + + valid_pix = mask > 0 + Npix = valid_pix.sum() + + inds = np.random.choice(Npix, size=N, replace=False) + + pixs = enmap.pixmap(shape, wcs) + coords = mask.pix2sky(pixs[:, valid_pix][:, inds], safe=False) / utils.degree + + # cmapper = catalogs.CatMapper(coords[1], coords[0], shape=shape, wcs=wcs) + # io.hplot( + # enmap.downgrade(cmapper.counts, 16), + # f"{paths.data}mask/randcounts_{choice}", + # mask=0, + # ) + + io.save_cols( + output_path + f"{sim_type}_{choice}_randoms.txt", + (coords[1], coords[0]), # ra, dec + ) + + + +def apply_snr_cut(catalog, threshold, label): + + keep = catalog["snrs"] > threshold + + zs = catalog["zs"][keep] + ras = catalog["ras"][keep] + decs = catalog["decs"][keep] + m200 = catalog["true_M200c"][keep] + m500 = catalog["true_M500c"][keep] + snrs = catalog["snrs"][keep] + + print(f"\nSNR > {threshold}") + print(f"N={len(zs)} [min SNR={snrs.min():.1f}]") + print(f"mean true M200c={m200.mean():.2f}") + print(f"mean true M500c={m500.mean():.2f}") + print(f"mean z={zs.mean():.2f}") + + io.save_cols( + output_path + f"{sim_type}_tsz_{label}_true.txt", + (ras, decs, zs, m200), + ) + + return zs, m200 + + +def resample_halo(tsz_zs, tsz_masses, snrcut): + + halo_file = "/home3/eunseong/cmbhalolensing/data/agora_halo_5e13cut.fits" + hdu = fits.open(halo_file) + data = hdu[1].data + + ras = data["RADeg"] + decs = data["decDeg"] + masses = data["true_M200c"] + M500c = data["true_M500c"] + zs = data["redshift"] + + print("\n ::: resampling halo sample!") + print(" ::: number of clusters in tsz cat =", len(tsz_masses)) + + num_bins = 30 + mbins = np.linspace(tsz_masses.min(), tsz_masses.max(), num_bins) + zbins = np.linspace(tsz_zs.min(), tsz_zs.max(), num_bins) + + hist, mbin_edges, zbin_edges = np.histogram2d( + tsz_masses, tsz_zs, bins=[mbins, zbins] + ) + + seed = 101 + rng = np.random.default_rng(seed) + + out = {k: [] for k in ["ras", "decs", "masses", "M500c", "zs"]} + + for i in range(len(mbin_edges) - 1): + for j in range(len(zbin_edges) - 1): + mmin, mmax = mbin_edges[i], mbin_edges[i + 1] + zmin, zmax = zbin_edges[j], zbin_edges[j + 1] + + mask = ( + (masses > mmin) + & (masses <= mmax) + & (zs > zmin) + & (zs <= zmax) + ) + + if hist[i, j] == 0 or not np.any(mask): + continue + + idx_pool = np.where(mask)[0] + n_target = int(hist[i, j]) + + replace = len(idx_pool) < n_target + inds = rng.choice(idx_pool, size=n_target, replace=replace) + + out["ras"].append(ras[inds]) + out["decs"].append(decs[inds]) + out["masses"].append(masses[inds]) + out["M500c"].append(M500c[inds]) + out["zs"].append(zs[inds]) + + ras = np.concatenate(out["ras"]) + decs = np.concatenate(out["decs"]) + masses = np.concatenate(out["masses"]) + M500c = np.concatenate(out["M500c"]) + zs = np.concatenate(out["zs"]) + + print(f" ::: z range = {zs.min():.2f}–{zs.max():.2f} (mean={zs.mean():.2f})") + print(f" ::: M200c range = {masses.min():.2f}–{masses.max():.2f} (mean={masses.mean():.2f})") + print(f" ::: M500c range = {M500c.min():.2f}–{M500c.max():.2f} (mean={M500c.mean():.2f})") + print(f" ::: number of halos after resampling = {len(ras)}") + + if snrcut == "snr4": + io.save_cols( + output_path + f"{sim_type}_halo_true.txt", + (ras, decs, zs, masses), + ) + else: + io.save_cols( + output_path + f"{sim_type}_halo_{snrcut}_true.txt", + (ras, decs, zs, masses), + ) + + + + +file_name = "data/20260320_cmb_tsz_ksz_mass_true.fits" # nsk61 +catalog = load_cat(file_name) + +generate_random_cat(n_per_obj=100, ras=catalog["ras"], choice="tsz") +generate_random_cat(n_per_obj=100, ras=catalog["ras"], choice="halo") + + +z_snr5p5, m_snr5p5 = apply_snr_cut(catalog, 5.5, "snr5p5") +z_snr7, m_snr7 = apply_snr_cut(catalog, 7.0, "snr7") + + +resample_halo(catalog["zs"], catalog["true_M200c"], "snr4") +resample_halo(z_snr5p5, m_snr5p5, "snr5p5") +resample_halo(z_snr7, m_snr7, "snr7") \ No newline at end of file diff --git a/input/sim_data.yml b/input/sim_data.yml index a7bc74b..e6c083d 100644 --- a/input/sim_data.yml +++ b/input/sim_data.yml @@ -7,32 +7,39 @@ sehgal_output_path: /home3/eunseong/cmbhalolensing/sim_results/sehgal agora_output_path: /home3/eunseong/cmbhalolensing/sim_results/agora simsuite_path: /home3/eunseong/cmbhalolensing/data/simsuite cat_path: /home3/eunseong/cmbhalolensing/data/cats +extra_map_path: /home3/eunseong/cmbhalolensing/data/lensed_sim_map # nemo sim catalogues -websky_tsz_cat: /home3/eunseong/nemo-sims/nemo-sim-kit11/sim-kit_NemoWebSky_CustomLensedCMB_DR5White_tenToA0Tuned_ACT-DR5-2Pass-GaussBeams_2degTiles/NemoWebSky_CustomLensedCMB_DR5White_tenToA0Tuned_ACT-DR5-2Pass-GaussBeams_2degTiles/NemoWebSky_CustomLensedCMB_DR5White_tenToA0Tuned_ACT-DR5-2Pass-GaussBeams_2degTiles_mass.fits -sehgal_tsz_cat: /home3/eunseong/nemo-sims/nemo-sim-kit13/sim-kit_NemoSehgal_CustomLensedCMB_DR5White_tenToA0Tuned_ACT-DR5-2Pass-GaussBeams_2degTiles/NemoSehgal_CustomLensedCMB_DR5White_tenToA0Tuned_ACT-DR5-2Pass-GaussBeams_2degTiles/NemoSehgal_CustomLensedCMB_DR5White_tenToA0Tuned_ACT-DR5-2Pass-GaussBeams_2degTiles_mass.fits -agora_tsz_cat: /home3/eunseong/nemo-sims/nemo-sim-kit11/sim-kit_NemoWebSky_CustomLensedCMB_DR5White_tenToA0Tuned_ACT-DR5-2Pass-GaussBeams_2degTiles/NemoWebSky_CustomLensedCMB_DR5White_tenToA0Tuned_ACT-DR5-2Pass-GaussBeams_2degTiles/NemoWebSky_CustomLensedCMB_DR5White_tenToA0Tuned_ACT-DR5-2Pass-GaussBeams_2degTiles_mass.fits +websky_tsz_cat: /home3/eunseong/nemo-sims/nemo-sim-kit54/sim-kit_websky_cmb_tsz_test1/websky_cmb_tsz_test1/websky_cmb_tsz_test1_mass.fits +sehgal_tsz_cat: /home3/eunseong/nemo-sims/nemo-sim-kit57/sehgal_cmb_tsz/sehgal_cmb_tsz_mass.fits +# agora_tsz_cat: /home3/eunseong/nemo-sims/nemo-sim-kit48/agora_cmb_tsz/agora_cmb_tsz_mass.fits +agora_tsz_cat: /home3/eunseong/nemo-sims/nemo-sim-kit61/20260320_cmb_tsz_ksz/20260320_cmb_tsz_ksz/20260320_cmb_tsz_ksz_mass.fits +# agora_tsz_cat: /home3/eunseong/nemo-sims/nemo-sim-kit62/20260320_cmb_tsz_cib/20260320_cmb_tsz_cib/20260320_cmb_tsz_cib_mass.fits +# agora_tsz_cat: /home3/eunseong/nemo-sims/nemo-sim-kit63/20260320_cmb_tsz_ksz_cib/20260320_cmb_tsz_ksz_cib/20260320_cmb_tsz_ksz_cib_mass.fits -nemosim_version: nemo-sim-kit11 -# nemo-sim-kit3: websky lensed_alm.fits + dr5white + 1/f + actbeam -# nemo-sim-kit4: websky lensed_alm.fits + flatwhite + actbeam -# nemo-sim-kit5: websky dlensed.fits + dr5white + 1/f + actbeam -# nemo-sim-kit7: websky dlensed.fits + flatwhite + actbeam -# nemo-sim-kit9: sehgal + flatwhite + actbeam -# (nemo updated 2024.07.03) -# nemo-sim-kit10: websky dlensed.fits + flatwhite + gaussbeam -# nemo-sim-kit11: websky dlensed.fits + dr5white + gaussbeam -# nemo-sim-kit12: sehgal + flatwhite + gaussbeam -# nemo-sim-kit13: sehgal + dr5white + gaussbeam +nemosim_version: nemo-sim-kit61 +# nemo-sim-kit47: agora from Matt (dr5white orphics + gaussbeam) +# nemo-sim-kit48: agora from Matt (dr5white orphics + gaussbeam) beam corrected - labeled as take2 from Matt +# nemo-sim-kit61: agora CMB+tSZ+kSZ (from Matt) 2026.03.23 +# nemo-sim-kit62: agora CMB+tSZ+CIB (from Matt) 2026.03.23 +# nemo-sim-kit63: agora CMB+tSZ+kSZ+CIB (from Matt) 2026.03.24 -simsuite_version: v0 -# websky_v0 (lensed_alm.fits) -# websky_v1 (dlensed.fits) -# sehgal_v0 -# agora_v0 +simsuite_version: v10 +# agora_v0 +# agora_v1: different noise seed used for ocmb_tsz +# agora_v2: planck beam and noise for grad (ocmb) +# agora_v3: planck beam and noise for grad (ocmb) and noise seeds are updated +# agora_v4: tsz map is switched from lensed (bahama78) to unlensed (bahama80) +# agora_v5: just organising and renaming - nothing changed! +# agora_v6: cmb only maps without noise and beam for the test +# agora_v7: resetting seeds for noise, checking with ksz and cib included +# agora_v8: for CMB+tSZ test +# agora_v9: copy of v4 +# so far all maps are full sky (later actfoot size maps are added for model subtraction) + # reprojected websky maps websky_kappa_reproj: reproj/kap_reproj.fits websky_cmb_reproj: reproj/lensed_alm_reproj.fits @@ -40,7 +47,9 @@ websky_tsz_reproj: reproj/tsz_8192_reproj.fits websky_kappa4p5_reproj: reproj/kap_lt4.5_reproj.fits websky_dlensed_reproj: reproj/dlensed_reproj.fits -websky_dtsz_reproj: reproj/dtsz_8192_reproj.fits +# websky_dtsz_reproj: reproj/dtsz_8192_reproj.fits +websky_dtsz090_reproj: reproj/dtsz090_8192_reproj.fits +websky_dtsz150_reproj: reproj/dtsz150_8192_reproj.fits websky_dksz_reproj: reproj/dksz_reproj.fits websky_dcib093_reproj: reproj/dcib_nu0093_reproj.fits websky_dcib145_reproj: reproj/dcib_nu0145_reproj.fits @@ -50,19 +59,19 @@ sehgal_kappa_reproj: reproj/healpix_4096_KappaeffLSStoCMBfullsky_reproj.fits sehgal_cmb_reproj: reproj/Sehgalsimparams_healpix_4096_KappaeffLSStoCMBfullsky_phi_SimLens_Tsynfastnopell_fast_lmax8000_nside4096_interp2.5_method1_1_lensed_map_reproj.fits sehgal_tsz_reproj: reproj/tSZ_skymap_healpix_nopell_Nside4096_y_tSZrescale0p75_reproj.fits -# reprojected agora maps -agora_kappa_reproj: ../../reproj/kappa_alm_recon_agora_phiNG_phi1_seed1_reproj.fits -agora_cmb_reproj: ../../reproj/agora_tqu1_phiNG_seed1_lmax16000_nside8192_interp1.6_method1_pol_1_lensedmap_reproj.fits -agora_tsz_reproj: ../../reproj/agora_ltszNG_bahamas78_bnd_unb_1.0e+12_1.0e+18_lensed_reproj.fits - +sehgal_dkappa_reproj: reproj/dhealpix_4096_KappaeffLSStoCMBfullsky_reproj.fits +sehgal_dlensed_reproj: reproj/dSehgalsimparams_healpix_4096_KappaeffLSStoCMBfullsky_phi_SimLens_Tsynfastnopell_fast_lmax8000_nside4096_interp2.5_method1_1_lensed_map_reproj.fits +# sehgal_dtsz_reproj: reproj/dtSZ_skymap_healpix_nopell_Nside4096_y_tSZrescale0p75_reproj.fits +sehgal_dtsz090_reproj: reproj/dtSZ090_skymap_healpix_nopell_Nside4096_y_tSZrescale0p75_reproj.fits +sehgal_dtsz150_reproj: reproj/dtSZ150_skymap_healpix_nopell_Nside4096_y_tSZrescale0p75_reproj.fits -# websky_tsz_cat: /home3/eunseong/nemo-sims/nemo-sim-kit7/sim-kit_NemoWebSky_CustomLensedCMB_tenToA0Tuned_ACT-DR5-2Pass_2degTiles/NemoWebSky_CustomLensedCMB_tenToA0Tuned_ACT-DR5-2Pass_2degTiles/NemoWebSky_CustomLensedCMB_tenToA0Tuned_ACT-DR5-2Pass_2degTiles_mass.fits -# sehgal_tsz_cat: /home3/eunseong/nemo-sims/nemo-sim-kit9/sim-kit_NemoSehgal_CustomLensedCMB_tenToA0Tuned_ACT-DR5-2Pass_2degTiles/NemoSehgal_CustomLensedCMB_tenToA0Tuned_ACT-DR5-2Pass_2degTiles/NemoSehgal_CustomLensedCMB_tenToA0Tuned_ACT-DR5-2Pass_2degTiles_mass.fits -# sehgal_tsz_cat: /home3/eunseong/nemo-sims/nemo-sim-kit12/sim-kit_NemoSehgal_CustomLensedCMB_flatWhite_tenToA0Tuned_ACT-DR5-2Pass-GaussBeams_2degTiles/NemoSehgal_CustomLensedCMB_flatWhite_tenToA0Tuned_ACT-DR5-2Pass-GaussBeams_2degTiles/NemoSehgal_CustomLensedCMB_flatWhite_tenToA0Tuned_ACT-DR5-2Pass-GaussBeams_2degTiles_mass.fits -# agora_tsz_cat: /home3/eunseong/nemo-sims/nemo-sim-kit10/sim-kit_NemoWebSky_CustomLensedCMB_flatWhite_tenToA0Tuned_ACT-DR5-2Pass-GaussBeams_2degTiles/NemoWebSky_CustomLensedCMB_flatWhite_tenToA0Tuned_ACT-DR5-2Pass-GaussBeams_2degTiles/NemoWebSky_CustomLensedCMB_flatWhite_tenToA0Tuned_ACT-DR5-2Pass-GaussBeams_2degTiles_mass.fits +# reprojected agora maps +agora_kappa_reproj: ../../es_to_delete/reproj/kappa_alm_recon_agora_phiNG_phi1_seed1_reproj.fits +agora_cmb_reproj: ../../reproj/agora_tqu1_phiNG_seed1_lmax16000_nside8192_interp1.6_method1_pol_1_lensedmap_reproj.fits +# agora_tsz_reproj: ../../reproj/agora_ltszNG_bahamas78_bnd_unb_1.0e+12_1.0e+18_lensed_reproj.fits +# agora_tsz_reproj: ../../reproj/agora_utszNG_bahamas80_bnd_unb_1.0e+12_1.0e+18_reproj.fits # unlensed tSZ +agora_tsz090_reproj: ../../reproj/agora_utszNG090_bahamas80_bnd_unb_1.0e+12_1.0e+18_reproj.fits # unlensed tSZ +agora_tsz150_reproj: ../../reproj/agora_utszNG150_bahamas80_bnd_unb_1.0e+12_1.0e+18_reproj.fits # unlensed tSZ -# # beams (corresponds to the beam used in NEMO) -# actbeam_f150: /data5/act/maps/dr5/beams/act_planck_dr5.01_s08s18_f150_daynight_beam.txt -# actbeam_f090: /data5/act/maps/dr5/beams/act_planck_dr5.01_s08s18_f090_daynight_beam.txt diff --git a/make_sim_suite_act.py b/make_sim_suite_act.py new file mode 100644 index 0000000..3ddc0a1 --- /dev/null +++ b/make_sim_suite_act.py @@ -0,0 +1,348 @@ +# sim_suite for model image subtracion run +# this should be done in ACT footprint maps (otherwise, cannot be subtracted) + +import time as t +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np +from pixell import enmap, reproject, utils, curvedsky, wcsutils, bunch +from orphics import maps, io +import healpy as hp +from past.utils import old_div +import argparse +import sys + +start = t.time() + +paths = bunch.Bunch(io.config_from_yaml("input/sim_data.yml")) +print("Paths: ", paths) + + +# RUN SETTING ------------------------------------------------------------------ + +parser = argparse.ArgumentParser() +parser.add_argument( + "save_name", type=str, help="Name you want for your output." +) +parser.add_argument( + "which_sim", type=str, help="Choose the sim e.g. websky or sehgal or agora." +) + +args = parser.parse_args() + + +output_path = paths.simsuite_path +save_name = args.save_name +if args.which_sim == "websky": + print(" ::: producing maps for WEBSKY sim") + save_name = "websky_" + save_name + sim_path = paths.websky_sim_path +elif args.which_sim == "sehgal": + print(" ::: producing maps for SEHGAL sim") + save_name = "sehgal_" + save_name + sim_path = paths.sehgal_sim_path +elif args.which_sim == "agora": + print(" ::: producing maps for AGORA sim") + save_name = "agora_" + save_name + sim_path = paths.agora_sim_path +save_dir = f"{output_path}/{save_name}/" +io.mkdir(f"{save_dir}") +print(" ::: saving to", save_dir) + + +# SIM SETTING ------------------------------------------------------------------ + +px = 0.5 + +fwhm_plc = 5.0 + +fwhm150 = 1.5 +fwhm090 = 2.2 + +nlevel = 15.0 +nlevel_plc = 35.0 + + +# BEAM CONVOLUTION ------------------------------------------------------------- + +def apply_beam(imap, ifwhm): + # map2alm of the maps, almxfl(alm, beam_1d) to convolve with beam, alm2map to convert back to map + if args.which_sim == "websky" or args.which_sim == "agora": nside = 8192 + elif args.which_sim == "sehgal": nside = 4096 + alm_lmax = nside * 3 + bfunc = lambda x: maps.gauss_beam(ifwhm, x) + imap_alm = curvedsky.map2alm(imap, lmax=alm_lmax) + beam_convoloved_alm = curvedsky.almxfl(imap_alm, bfunc) + return curvedsky.alm2map(beam_convoloved_alm, enmap.empty(imap.shape, imap.wcs)) + +# PREPARING MAPS --------------------------------------------------------------- + +# reading lensed cmb map + +ifile = f"{paths.mat_path}/dlensed_actfoot.fits" # act footprint +lmap = enmap.read_map(ifile) +print(" ::: reading a map for ACT footprint geometry:", ifile) +shape, wcs = lmap.shape, lmap.wcs + + +if args.which_sim == "agora": + ifile = f"{sim_path}cmb/len/tqu1/agora_tqu1_phiNG_seed1_lmax16000_nside8192_interp1.6_method1_pol_1_lensedmap.fits" + hmap,_,_ = hp.read_map(ifile, field=[0,1,2]).astype(np.float64) + lmap = reproject.healpix2map(hmap, shape, wcs) + print(" ::: reading lensed cmb map:", ifile) + +print("lmap", lmap.shape) +shape, wcs = lmap.shape, lmap.wcs + + + +# reading tsz map + +if args.which_sim == "agora": + ifile = f"{sim_path}/tsz/unl/agora_utszNG_bahamas80_bnd_unb_1.0e+12_1.0e+18.fits" + hmap = hp.read_map(ifile).astype(np.float64) + ymap = reproject.healpix2map(hmap, shape, wcs) + +print(" ::: reading ymap:", ifile) +print("ymap", ymap.shape) + +# convert compton-y to delta-T (in uK) +tcmb = 2.726 +tcmb_uK = tcmb * 1e6 #micro-Kelvin +H_cgs = 6.62608e-27 +K_cgs = 1.3806488e-16 + +def fnu(nu): + """ + nu in GHz + tcmb in Kelvin + """ + mu = H_cgs*(1e9*nu)/(K_cgs*tcmb) + ans = mu/np.tanh(old_div(mu,2.0)) - 4.0 + return ans + +freq_sz = 90 +print(" ::: map frequency at %d GHz" %freq_sz) +tszmap090 = fnu(freq_sz) * ymap * tcmb_uK +print(" ::: converting ymap to tsz map at %d GHz" %freq_sz) +print("tszmap090", tszmap090.shape) + +freq_sz = 150 +print(" ::: map frequency at %d GHz" %freq_sz) +tszmap150 = fnu(freq_sz) * ymap * tcmb_uK +print(" ::: converting ymap to tsz map at %d GHz" %freq_sz) +print("tszmap150", tszmap150.shape) + + +# reading ksz map + +# r_kszmap = sim_path + paths.agora_ksz_reproj +if args.which_sim == "agora": + ifile = f"{sim_path}ksz/unl/agora_ukszNG_bahamas80_bnd_unb_1.0e+12_1.0e+18.fits" + hmap = hp.read_map(ifile).astype(np.float64) + kszmap = reproject.healpix2map(hmap, shape, wcs) + +# enmap.write_map(r_kszmap, kszmap) +print(" ::: reading and reprojecting ksz map:", ifile) +print("kszmap", kszmap.shape) + + + +# reading cib map +if args.which_sim == "agora": + + fluxcut_val = 5. # flux per pixel [mJy] + + def flux_density_to_temp(freq_GHz): + # get factor for converting delta flux density in [MJy/sr] to delta T in CMB units [uK] + freq = float(freq_GHz) + x = freq / 56.8 + return (1.05e3 * (np.exp(x)-1)**2 * + np.exp(-x) * (freq / 100)**-4) + + freq = 150 + ifile = f"{sim_path}cib/uK/len/act/agora_len_mag_cibmap_act_{freq}ghz_uk.fits" # in uK + hmap = hp.read_map(ifile).astype(np.float64) + print(hmap.min(), hmap.max(), hmap.mean()) + + nside = hp.get_nside(hmap) # extract the HEALPix resolution + print(nside) + + pixel_solid_angle = hp.nside2pixarea(nside) # in steradian + lim = fluxcut_val * 1e-9 / pixel_solid_angle # [mJy] to [MJy/sr] + uK_lim = flux_density_to_temp(freq) * lim # [MJy/sr] to [uK] + + print(f"conversion factor: {flux_density_to_temp(freq):.2f}") + print(f"uK threshold: {uK_lim:.2f}") + print(f"number of pixels above threshold: {np.sum(hmap > uK_lim)} ({np.sum(hmap > uK_lim)/np.sum(hmap)*100:.6f} percent)") + print(f"implementing flux cut on map of {fluxcut_val} mJy per pixel") + hmap[hmap > uK_lim]=0. + print(hmap.min(), hmap.max(), hmap.mean()) + + cibmap150 = reproject.healpix2map(hmap, shape, wcs, method="spline") + + + freq = 90 + ifile = f"{sim_path}cib/uK/len/act/agora_len_mag_cibmap_act_{freq}ghz_uk.fits" # in uK + hmap = hp.read_map(ifile).astype(np.float64) + print(hmap.min(), hmap.max(), hmap.mean()) + + nside = hp.get_nside(hmap) # extract the HEALPix resolution + print(nside) + + pixel_solid_angle = hp.nside2pixarea(nside) # in steradian + lim = fluxcut_val * 1e-9 / pixel_solid_angle # [mJy] to [MJy/sr] + uK_lim = flux_density_to_temp(freq) * lim # [MJy/sr] to [uK] + + print(f"conversion factor: {flux_density_to_temp(freq):.2f}") + print(f"uK threshold: {uK_lim:.2f}") + print(f"number of pixels above threshold: {np.sum(hmap > uK_lim)} ({np.sum(hmap > uK_lim)/np.sum(hmap)*100:.6f} percent)") + print(f"implementing flux cut on map of {fluxcut_val} mJy per pixel") + hmap[hmap > uK_lim]=0. + print(hmap.min(), hmap.max(), hmap.mean()) + + + cibmap090 = reproject.healpix2map(hmap, shape, wcs, method="spline") + + # enmap.write_map(r_cibmap, cibmap) + print("cibmap090", cibmap090.shape) + print("cibmap150", cibmap150.shape) + + +# assert wcsutils.equal(lmap.wcs, kmap.wcs) +assert wcsutils.equal(lmap.wcs, tszmap090.wcs) +assert wcsutils.equal(lmap.wcs, tszmap150.wcs) +assert wcsutils.equal(lmap.wcs, kszmap.wcs) +assert wcsutils.equal(lmap.wcs, cibmap090.wcs) +assert wcsutils.equal(lmap.wcs, cibmap150.wcs) + + + + + +scmb_tsz150 = lmap + tszmap150 +scmb_tsz090 = lmap + tszmap090 +scmb_tsz150_ksz = lmap + tszmap150 + kszmap +scmb_tsz090_ksz = lmap + tszmap090 + kszmap +scmb_tsz150_cib = lmap + tszmap150 + cibmap150 +scmb_tsz090_cib = lmap + tszmap090 + cibmap090 +scmb_tsz150_ksz_cib = lmap + tszmap150 + kszmap + cibmap150 +scmb_tsz090_ksz_cib = lmap + tszmap090 + kszmap + cibmap090 + +print(" ::: SIGNAL maps are ready!") +print(" ::: applying beam and adding white noise") + +h_ocmb_tsz150 = apply_beam(scmb_tsz150, fwhm150) + maps.white_noise(shape, wcs, noise_muK_arcmin=nlevel, seed=12) +h_ocmb_tsz090 = apply_beam(scmb_tsz090, fwhm090) + maps.white_noise(shape, wcs, noise_muK_arcmin=nlevel, seed=13) +h_ocmb_tsz150_ksz = apply_beam(scmb_tsz150_ksz, fwhm150) + maps.white_noise(shape, wcs, noise_muK_arcmin=nlevel, seed=12) +h_ocmb_tsz090_ksz = apply_beam(scmb_tsz090_ksz, fwhm090) + maps.white_noise(shape, wcs, noise_muK_arcmin=nlevel, seed=13) +h_ocmb_tsz150_cib = apply_beam(scmb_tsz150_cib, fwhm150) + maps.white_noise(shape, wcs, noise_muK_arcmin=nlevel, seed=12) +h_ocmb_tsz090_cib = apply_beam(scmb_tsz090_cib, fwhm090) + maps.white_noise(shape, wcs, noise_muK_arcmin=nlevel, seed=13) +h_ocmb_tsz150_ksz_cib = apply_beam(scmb_tsz150_ksz_cib, fwhm150) + maps.white_noise(shape, wcs, noise_muK_arcmin=nlevel, seed=12) +h_ocmb_tsz090_ksz_cib = apply_beam(scmb_tsz090_ksz_cib, fwhm090) + maps.white_noise(shape, wcs, noise_muK_arcmin=nlevel, seed=13) + + + +print(" ::: OBSERVED maps are ready!") + +# save_h_ocmb_tsz150_ksz = f"{save_dir}h_ocmb_tsz150_ksz_actfoot.fits" +# save_h_ocmb_tsz090_ksz = f"{save_dir}h_ocmb_tsz090_ksz_actfoot.fits" +# enmap.write_map(save_h_ocmb_tsz150_ksz, h_ocmb_tsz150_ksz) +# enmap.write_map(save_h_ocmb_tsz090_ksz, h_ocmb_tsz090_ksz) +# sys.exit() + + + + +# reading the model map -- model map is beam convolved + +## only tSZ +# model150 = f"/home3/eunseong/nemo-sims/nemo-sim-kit48/agora_cmb_tsz/clusterModelMaps/clusterModelMap_f150.fits" +# model090 = f"/home3/eunseong/nemo-sims/nemo-sim-kit48/agora_cmb_tsz/clusterModelMaps/clusterModelMap_f090.fits" +## tSZ with kSZ +model150 = f"/home3/eunseong/nemo-sims/nemo-sim-kit61/20260320_cmb_tsz_ksz/20260320_cmb_tsz_ksz/clusterModelMaps/clusterModelMap_f150.fits" +model090 = f"/home3/eunseong/nemo-sims/nemo-sim-kit61/20260320_cmb_tsz_ksz/20260320_cmb_tsz_ksz/clusterModelMaps/clusterModelMap_f090.fits" +## tSZ with CIB +# model150 = f"/home3/eunseong/nemo-sims/nemo-sim-kit62/20260320_cmb_tsz_cib/20260320_cmb_tsz_cib/clusterModelMaps/clusterModelMap_f150.fits" +# model090 = f"/home3/eunseong/nemo-sims/nemo-sim-kit62/20260320_cmb_tsz_cib/20260320_cmb_tsz_cib/clusterModelMaps/clusterModelMap_f090.fits" +## tSZ with kSZ+CIB +# model150 = f"/home3/eunseong/nemo-sims/nemo-sim-kit63/20260320_cmb_tsz_ksz_cib/20260320_cmb_tsz_ksz_cib/clusterModelMaps/clusterModelMap_f150.fits" +# model090 = f"/home3/eunseong/nemo-sims/nemo-sim-kit63/20260320_cmb_tsz_ksz_cib/20260320_cmb_tsz_ksz_cib/clusterModelMaps/clusterModelMap_f090.fits" + + + +print(" ::: reading model image map:", model150) +print(" ::: reading model image map:", model090) + +modelmap150 = enmap.read_map(model150, delayed=False) +modelmap090 = enmap.read_map(model090, delayed=False) +print(modelmap150.shape, modelmap090.shape) + + + +print(" ::: subtracting the model maps!") + + +# h_ocmb_msub150 = h_ocmb_tsz150 - modelmap150 +# h_ocmb_msub090 = h_ocmb_tsz090 - modelmap090 +h_ocmb_msub150_ksz = h_ocmb_tsz150_ksz - modelmap150 +h_ocmb_msub090_ksz = h_ocmb_tsz090_ksz - modelmap090 +# h_ocmb_msub150_cib = h_ocmb_tsz150_cib - modelmap150 +# h_ocmb_msub090_cib = h_ocmb_tsz090_cib - modelmap090 +# h_ocmb_msub150_ksz_cib = h_ocmb_tsz150_ksz_cib - modelmap150 +# h_ocmb_msub090_ksz_cib = h_ocmb_tsz090_ksz_cib - modelmap090 + + +print(" ::: model subtracted maps are ready!") + + + +# save_h_ocmb_msub150 = f"{save_dir}h_ocmb_msub150_test.fits" +# save_h_ocmb_msub090 = f"{save_dir}h_ocmb_msub090_test.fits" +save_h_ocmb_msub150_ksz = f"{save_dir}h_ocmb_msub150_ksz_updated.fits" +save_h_ocmb_msub090_ksz = f"{save_dir}h_ocmb_msub090_ksz_updated.fits" +# save_h_ocmb_msub150_cib = f"{save_dir}h_ocmb_msub150_cib5_updated.fits" +# save_h_ocmb_msub090_cib = f"{save_dir}h_ocmb_msub090_cib5_updated.fits" +# save_h_ocmb_msub150_ksz_cib = f"{save_dir}h_ocmb_msub150_ksz_cib5_updated.fits" +# save_h_ocmb_msub090_ksz_cib = f"{save_dir}h_ocmb_msub090_ksz_cib5_updated.fits" + + + +# enmap.write_map(save_h_ocmb_msub150, h_ocmb_msub150) +# enmap.write_map(save_h_ocmb_msub090, h_ocmb_msub090) +enmap.write_map(save_h_ocmb_msub150_ksz, h_ocmb_msub150_ksz) +enmap.write_map(save_h_ocmb_msub090_ksz, h_ocmb_msub090_ksz) +# enmap.write_map(save_h_ocmb_msub150_cib, h_ocmb_msub150_cib) +# enmap.write_map(save_h_ocmb_msub090_cib, h_ocmb_msub090_cib) +# enmap.write_map(save_h_ocmb_msub150_ksz_cib, h_ocmb_msub150_ksz_cib) +# enmap.write_map(save_h_ocmb_msub090_ksz_cib, h_ocmb_msub090_ksz_cib) + + + +print(" ::: all maps are saved! yayyyyyy") + + + + + + + +elapsed = t.time() - start +print("\r ::: entire run took %.1f seconds" %elapsed) + + + + + + + + + + + + + + + + diff --git a/recon_sim.py b/recon_sim.py index cf56616..1a704a2 100644 --- a/recon_sim.py +++ b/recon_sim.py @@ -386,10 +386,12 @@ def bin(data, modrmap, bin_edges): if args.cmb_ksz_msub: - cmb_msub090 = f"{simsuite_path}h_ocmb_msub090_ksz_updated.fits" - cmb_msub150 = f"{simsuite_path}h_ocmb_msub150_ksz_updated.fits" + # cmb_msub090 = f"{simsuite_path}h_ocmb_msub090_ksz_updated.fits" + # cmb_msub150 = f"{simsuite_path}h_ocmb_msub150_ksz_updated.fits" + cmb_msub090 = f"{simsuite_path}/../../maps/agora_h_ocmb_tsz_ksz_f090_msub.fits" # new model subtraction method + cmb_msub150 = f"{simsuite_path}/../../maps/agora_h_ocmb_tsz_ksz_f150_msub.fits" g_cmb_ksz = f"{simsuite_path}g_ocmb_ksz.fits" - + print(" ::: reading the model subtracted maps!", cmb_msub090) print(" ::: reading the model subtracted maps!", cmb_msub150) @@ -1559,4 +1561,4 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): elapsed = t.time() - start -print("\r ::: entire run took %.1f seconds" %elapsed) +print("\r ::: entire run took %.1f seconds" %elapsed) \ No newline at end of file diff --git a/stack.py b/stack.py index fc5cf08..1fe53f3 100644 --- a/stack.py +++ b/stack.py @@ -9,7 +9,8 @@ from scipy.optimize import curve_fit from numpy import save import time -import symlens +# import symlens +from symlens.qe import QE, N_l, N_l_from_A_l_optimal import healpy as hp import os, sys from enlib import bench @@ -103,7 +104,8 @@ else: if not (args.ilc_maps): act_map = (paths.act_data + data_choice.hres_150) - famap_150 = enmap.read_map(act_map, delayed=False, sel=np.s_[0, ...]) + fmap = enmap.read_map(act_map, delayed=False) + famap_150 = fmap[0] if fmap.ndim == 3 else fmap else: act_map = (paths.act_data + data_choice.hres) famap_150 = enmap.read_map(act_map, delayed=False) @@ -125,7 +127,8 @@ amap_90 = enmap.read_map(f'{paths.fullsim_path}/af090_sim_{args.full_sim_index:06d}.fits') else: act_map = (paths.act_data + data_choice.hres_090) - famap_90 = enmap.read_map(act_map, delayed=False, sel=np.s_[0, ...]) + fmap = enmap.read_map(act_map, delayed=False) + famap_90 = fmap[0] if fmap.ndim == 3 else fmap # SZ cluster model image subtraction for 90 GHz if args.hres_grad: @@ -148,23 +151,20 @@ act_map = ( paths.act_data + data_choice.hres_150_night ) - namap_150 = enmap.read_map(act_map, delayed=False, sel=np.s_[0, ...]) + fmap = enmap.read_map(act_map, delayed=False) + namap_150 = fmap[0] if fmap.ndim == 3 else fmap act_map = ( paths.act_data + data_choice.hres_090_night ) - namap_90 = enmap.read_map(act_map, delayed=False, sel=np.s_[0, ...]) - + fmap = enmap.read_map(act_map, delayed=False) + namap_90 = fmap[0] if fmap.ndim == 3 else fmap null_map_150 = famap_150 - namap_150 null_map_90 = famap_90 - namap_90 # Inv var map for 90 GHz ivar_map = (paths.act_data + data_choice.hres_ivar) - - try: - imap_90 = enmap.read_map(ivar_map, delayed=False, sel=np.s_[0, ...]) - except: - imap_90 = enmap.read_map(ivar_map, delayed=False) - + fmap = enmap.read_map(ivar_map, delayed=False) + imap_90 = fmap[0] if fmap.ndim == 3 else fmap rms_map = maps.rms_from_ivar( imap_90, cylindrical=True, safe=False # safe=False for DR6 maps @@ -1126,7 +1126,7 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): assert np.all(np.isfinite(feed_dict[key])) # ask for reconstruction in Fourier space - cqe = symlens.QE( + cqe = QE( shape, wcs, feed_dict, @@ -1236,12 +1236,12 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): if args.full_nl: # FIXME: The full calculation seems to sometimes be negative - Nl = symlens.N_l(shape,wcs,feed_dict,estimator="hdv_curl" if args.curl else "hdv",XY='TT', + Nl = N_l(shape,wcs,feed_dict,estimator="hdv_curl" if args.curl else "hdv",XY='TT', xmask=xmask,ymask=ymask, Al=Al,field_names=['P','A'],kmask=kmask,power_name="t") # !!! else: # Use approximate noise assuming estimator is optimal (Narrator: it isn't) - Nl = symlens.N_l_from_A_l_optimal(shape, wcs, Al) + Nl = N_l_from_A_l_optimal(shape, wcs, Al) cents, bnl = lbinner.bin(Nl) nmean = bnl[np.logical_and(cents > defaults.kappa_noise_mean_Lmin, cents < defaults.kappa_noise_mean_Lmax)].mean() if args.debug_powers: @@ -1498,4 +1498,3 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): elapsed = time.time() - start_time print("\r ::: entire run took %.1f seconds" % elapsed) - diff --git a/tsz_subtract.py b/tsz_subtract.py index 980b25a..41db01c 100644 --- a/tsz_subtract.py +++ b/tsz_subtract.py @@ -200,7 +200,8 @@ def deltaT_profile(m500c, z, cosmo, nu_ghz, ntheta=256, los_max_r500=5.0, los_np pe = p500 * gnfw_pressure(x) / P_TH_TO_PE # electron pressure [keV cm^-3] # y = (sigma_T / m_e c^2) * 2 * int_0^Lmax P_e dl, with dl in cm. prefac = SIGMA_T_CM2 / ME_C2_KEV - y_theta = prefac * 2.0 * np.trapezoid(pe, los * MPC_CM, axis=1) + # y_theta = prefac * 2.0 * np.trapezoid(pe, los * MPC_CM, axis=1) + y_theta = prefac * 2.0 * np.trapz(pe, los * MPC_CM, axis=1) # numpy version issue g_nu = foregrounds.g_tsz(nu_ghz) dT = g_nu * foregrounds.TCMB_uK * y_theta # micro-Kelvin @@ -236,7 +237,8 @@ def bfn(ell): return bfn except (TypeError, ValueError): - ells, bells = np.loadtxt(beam_arg, unpack=True) + # ells, bells = np.loadtxt(beam_arg, unpack=True) + ells, bells = np.loadtxt(beam_arg, usecols=[0,1], unpack=True) bells = bells / bells[np.argmin(ells)] interp = interp1d(ells, bells, bounds_error=False, fill_value=(bells[0], 0.0)) @@ -1058,7 +1060,7 @@ def log(msg): "ivar geometries differ." ) return - + # --- Subtract into full map and write ------------------------------------ log("Loading full coadd for subtraction...") t_load = time.time() @@ -1674,4 +1676,4 @@ def main(): if __name__ == "__main__": - main() + main() \ No newline at end of file From 92bd9f9a91d940b189545961642a51ddb96f55fa Mon Sep 17 00:00:00 2001 From: Eunseong Lee <36232273+eunseongleee@users.noreply.github.com> Date: Tue, 7 Jul 2026 13:34:04 -0400 Subject: [PATCH 14/14] another check point --- input/data.yml | 75 +++++-- input/sim_data.yml | 97 +++----- make_sim_suite.py | 542 +++++++++++++++++++++------------------------ recon_sim.py | 163 ++++++-------- 4 files changed, 411 insertions(+), 466 deletions(-) diff --git a/input/data.yml b/input/data.yml index 19fc961..5093854 100644 --- a/input/data.yml +++ b/input/data.yml @@ -1,28 +1,63 @@ # choose coadd maps to use for hres leg -hres_090: 20211219/act_planck_daynight_f090_map_srcfree.fits -hres_150: 20211219/act_planck_daynight_f150_map_srcfree.fits -hres_ivar: 20211219/act_planck_daynight_f090_ivar.fits -hres_model_090: 20230721/modelImages/DR6_candidates_SNR5_model_f090.fits -hres_model_150: 20230721/modelImages/DR6_candidates_SNR5_model_f150.fits +# # dr5 run! +# hres_090: maps/dr5/act_planck_dr5.01_s08s18_AA_f090_night_map_srcfree.fits +# hres_150: maps/dr5/act_planck_dr5.01_s08s18_AA_f150_night_map_srcfree.fits +# hres_ivar: maps/dr5/act_planck_dr5.01_s08s18_AA_f090_night_ivar.fits +# hres_model_090: hiltonm/models/dr5/S18d_202006_pass2_confirmed_model_f090.fits +# hres_model_150: hiltonm/models/dr5/S18d_202006_pass2_confirmed_model_f150.fits + + +# dr6 run! +# hres_090: coadds/20240323_simple/act_planck_s08_s22_f090_daynight_map_srcfree.fits +# hres_150: coadds/20240323_simple/act_planck_s08_s22_f150_daynight_map_srcfree.fits +hres_ivar: coadds/20240323_simple/act_planck_s08_s22_f090_daynight_ivar.fits +# hres_model_090: hiltonm/cluster_candidates/20230721/modelImages/DR6_candidates_SNR5_model_f090.fits +# hres_model_150: hiltonm/cluster_candidates/20230721/modelImages/DR6_candidates_SNR5_model_f150.fits +# hres_model_090: iabril/nemo_20250924/candidatesModelMap_f090_SNR5.fits +# hres_model_150: iabril/nemo_20250924/candidatesModelMap_f150_SNR5.fits +hres_090: ../../home3/eunseong/cmbhalolensing/data/maps/act_coadd_f090_msub.fits +hres_150: ../../home3/eunseong/cmbhalolensing/data/maps/act_coadd_f150_msub.fits +hres_090_night: coadds/20240323_simple/act_planck_s08_s22_f090_night_map_srcfree.fits +hres_150_night: coadds/20240323_simple/act_planck_s08_s22_f150_night_map_srcfree.fits + + +# # test ILC run! +# hres: ilc_dr6v3/20230606/ilc_fullRes_TT.fits +# grad: ilc_dr6v3/20230606/ilc_fullRes_deproj_sz_TT.fits +# hres_090: coadds/20211219/act_planck_daynight_f090_map_srcfree.fits +# hres_ivar: coadds/20211219/act_planck_daynight_f090_ivar.fits + +# # dr6 run! (this turns out to be old) +# hres_090: coadds/20220909_simple/act_planck_s08_s21_daynight_f090_map_srcfree.fits +# hres_150: coadds/20220909_simple/act_planck_s08_s21_daynight_f150_map_srcfree.fits +# hres_ivar: coadds/20220909_simple/act_planck_s08_s21_daynight_f090_ivar.fits +# hres_model_090: hiltonm/cluster_candidates/20230721/modelImages/DR6_candidates_SNR5_model_f090.fits +# hres_model_150: hiltonm/cluster_candidates/20230721/modelImages/DR6_candidates_SNR5_model_f150.fits + + + +# dr5_090: maps/dr5/act_planck_dr5.01_s08s18_AA_f090_night_map_srcfree.fits +# dr5_150: maps/dr5/act_planck_dr5.01_s08s18_AA_f150_night_map_srcfree.fits +# dr5_ivar: maps/dr5/act_planck_dr5.01_s08s18_AA_f090_night_ivar.fits + +# dr6_090: coadds/20211219/act_planck_daynight_f090_map_srcfree.fits +# dr6_150: coadds/20211219/act_planck_daynight_f150_map_srcfree.fits +# dr6_ivar: coadds/20211219/act_planck_daynight_f090_ivar.fits + +# dr6_simple_090: coadds/20220909_simple/act_planck_s08_s21_daynight_f090_map_srcfree.fits +# dr6_simple_150: coadds/20220909_simple/act_planck_s08_s21_daynight_f150_map_srcfree.fits +# dr6_simple_ivar: coadds/20220909_simple/act_planck_s08_s21_daynight_f090_ivar.fits + +# dr6_ilc: ilc_dr6v3/20230606/ilc_fullRes_TT.fits +# dr6_ilc_noSZ: ilc_dr6v3/20230606/ilc_fullRes_deproj_sz_TT.fits -# dr5_090: act_planck_dr5.01_s08s18_AA_f090_night_map_srcfree.fits -# dr5_150: act_planck_dr5.01_s08s18_AA_f150_night_map_srcfree.fits -# dr5_ivar: act_planck_dr5.01_s08s18_AA_f090_night_ivar.fits - -# dr6_090: 20211219/act_planck_daynight_f090_map_srcfree.fits -# dr6_150: 20211219/act_planck_daynight_f150_map_srcfree.fits -# dr6_ivar: 20211219/act_planck_daynight_f090_ivar.fits +# dr5_model_090: hiltonm/models/dr5/S18d_202006_pass2_confirmed_model_f090.fits +# dr5_model_150: hiltonm/models/dr5/S18d_202006_pass2_confirmed_model_f150.fits -# dr6_simple_090: 20220909_simple/act_planck_s08_s21_daynight_f090_map_srcfree.fits -# dr6_simple_150: 20220909_simple/act_planck_s08_s21_daynight_f150_map_srcfree.fits -# dr6_simple_ivar: 20220909_simple/act_planck_s08_s21_daynight_f090_ivar.fits - -# dr5_model_090: S18d_202006_pass2_confirmed_model_f090.fits -# dr5_model_150: S18d_202006_pass2_confirmed_model_f150.fits +# dr6_model_simple_090: hiltonm/cluster_candidates/20230721/modelImages/DR6_candidates_SNR5_model_f090.fits +# dr6_model_simple_150: hiltonm/cluster_candidates/20230721/modelImages/DR6_candidates_SNR5_model_f150.fits -# dr6_model_simple_090: 20230721/modelImages/DR6_candidates_SNR5_model_f090.fits -# dr6_model_simple_150: 20230721/modelImages/DR6_candidates_SNR5_model_f150.fits diff --git a/input/sim_data.yml b/input/sim_data.yml index e6c083d..c2aa524 100644 --- a/input/sim_data.yml +++ b/input/sim_data.yml @@ -7,71 +7,42 @@ sehgal_output_path: /home3/eunseong/cmbhalolensing/sim_results/sehgal agora_output_path: /home3/eunseong/cmbhalolensing/sim_results/agora simsuite_path: /home3/eunseong/cmbhalolensing/data/simsuite cat_path: /home3/eunseong/cmbhalolensing/data/cats +nsk_path: /home3/eunseong/nemo-sims/ extra_map_path: /home3/eunseong/cmbhalolensing/data/lensed_sim_map -# nemo sim catalogues -websky_tsz_cat: /home3/eunseong/nemo-sims/nemo-sim-kit54/sim-kit_websky_cmb_tsz_test1/websky_cmb_tsz_test1/websky_cmb_tsz_test1_mass.fits -sehgal_tsz_cat: /home3/eunseong/nemo-sims/nemo-sim-kit57/sehgal_cmb_tsz/sehgal_cmb_tsz_mass.fits -# agora_tsz_cat: /home3/eunseong/nemo-sims/nemo-sim-kit48/agora_cmb_tsz/agora_cmb_tsz_mass.fits -agora_tsz_cat: /home3/eunseong/nemo-sims/nemo-sim-kit61/20260320_cmb_tsz_ksz/20260320_cmb_tsz_ksz/20260320_cmb_tsz_ksz_mass.fits -# agora_tsz_cat: /home3/eunseong/nemo-sims/nemo-sim-kit62/20260320_cmb_tsz_cib/20260320_cmb_tsz_cib/20260320_cmb_tsz_cib_mass.fits -# agora_tsz_cat: /home3/eunseong/nemo-sims/nemo-sim-kit63/20260320_cmb_tsz_ksz_cib/20260320_cmb_tsz_ksz_cib/20260320_cmb_tsz_ksz_cib_mass.fits - -nemosim_version: nemo-sim-kit61 -# nemo-sim-kit47: agora from Matt (dr5white orphics + gaussbeam) -# nemo-sim-kit48: agora from Matt (dr5white orphics + gaussbeam) beam corrected - labeled as take2 from Matt -# nemo-sim-kit61: agora CMB+tSZ+kSZ (from Matt) 2026.03.23 -# nemo-sim-kit62: agora CMB+tSZ+CIB (from Matt) 2026.03.23 -# nemo-sim-kit63: agora CMB+tSZ+kSZ+CIB (from Matt) 2026.03.24 - - - -simsuite_version: v10 -# agora_v0 -# agora_v1: different noise seed used for ocmb_tsz -# agora_v2: planck beam and noise for grad (ocmb) -# agora_v3: planck beam and noise for grad (ocmb) and noise seeds are updated -# agora_v4: tsz map is switched from lensed (bahama78) to unlensed (bahama80) -# agora_v5: just organising and renaming - nothing changed! -# agora_v6: cmb only maps without noise and beam for the test -# agora_v7: resetting seeds for noise, checking with ksz and cib included -# agora_v8: for CMB+tSZ test -# agora_v9: copy of v4 -# so far all maps are full sky (later actfoot size maps are added for model subtraction) - -# reprojected websky maps -websky_kappa_reproj: reproj/kap_reproj.fits -websky_cmb_reproj: reproj/lensed_alm_reproj.fits -websky_tsz_reproj: reproj/tsz_8192_reproj.fits - -websky_kappa4p5_reproj: reproj/kap_lt4.5_reproj.fits -websky_dlensed_reproj: reproj/dlensed_reproj.fits -# websky_dtsz_reproj: reproj/dtsz_8192_reproj.fits -websky_dtsz090_reproj: reproj/dtsz090_8192_reproj.fits -websky_dtsz150_reproj: reproj/dtsz150_8192_reproj.fits -websky_dksz_reproj: reproj/dksz_reproj.fits -websky_dcib093_reproj: reproj/dcib_nu0093_reproj.fits -websky_dcib145_reproj: reproj/dcib_nu0145_reproj.fits - -# reprojected sehgal maps -sehgal_kappa_reproj: reproj/healpix_4096_KappaeffLSStoCMBfullsky_reproj.fits -sehgal_cmb_reproj: reproj/Sehgalsimparams_healpix_4096_KappaeffLSStoCMBfullsky_phi_SimLens_Tsynfastnopell_fast_lmax8000_nside4096_interp2.5_method1_1_lensed_map_reproj.fits -sehgal_tsz_reproj: reproj/tSZ_skymap_healpix_nopell_Nside4096_y_tSZrescale0p75_reproj.fits - -sehgal_dkappa_reproj: reproj/dhealpix_4096_KappaeffLSStoCMBfullsky_reproj.fits -sehgal_dlensed_reproj: reproj/dSehgalsimparams_healpix_4096_KappaeffLSStoCMBfullsky_phi_SimLens_Tsynfastnopell_fast_lmax8000_nside4096_interp2.5_method1_1_lensed_map_reproj.fits -# sehgal_dtsz_reproj: reproj/dtSZ_skymap_healpix_nopell_Nside4096_y_tSZrescale0p75_reproj.fits -sehgal_dtsz090_reproj: reproj/dtSZ090_skymap_healpix_nopell_Nside4096_y_tSZrescale0p75_reproj.fits -sehgal_dtsz150_reproj: reproj/dtSZ150_skymap_healpix_nopell_Nside4096_y_tSZrescale0p75_reproj.fits - - -# reprojected agora maps -agora_kappa_reproj: ../../es_to_delete/reproj/kappa_alm_recon_agora_phiNG_phi1_seed1_reproj.fits -agora_cmb_reproj: ../../reproj/agora_tqu1_phiNG_seed1_lmax16000_nside8192_interp1.6_method1_pol_1_lensedmap_reproj.fits -# agora_tsz_reproj: ../../reproj/agora_ltszNG_bahamas78_bnd_unb_1.0e+12_1.0e+18_lensed_reproj.fits -# agora_tsz_reproj: ../../reproj/agora_utszNG_bahamas80_bnd_unb_1.0e+12_1.0e+18_reproj.fits # unlensed tSZ -agora_tsz090_reproj: ../../reproj/agora_utszNG090_bahamas80_bnd_unb_1.0e+12_1.0e+18_reproj.fits # unlensed tSZ -agora_tsz150_reproj: ../../reproj/agora_utszNG150_bahamas80_bnd_unb_1.0e+12_1.0e+18_reproj.fits # unlensed tSZ +# agora_tsz_cat: agora_cmb_tsz/agora_cmb_tsz_mass.fits +agora_tsz_cat: 20260320_cmb_tsz_ksz/20260320_cmb_tsz_ksz/20260320_cmb_tsz_ksz_mass.fits +# agora_tsz_cat: 20260320_cmb_tsz_cib/20260320_cmb_tsz_cib/20260320_cmb_tsz_cib_mass.fits +# agora_tsz_cat: 20260320_cmb_tsz_ksz_cib/20260320_cmb_tsz_ksz_cib/20260320_cmb_tsz_ksz_cib_mass.fits +nemosim_version: nemo-sim-kit61 +# nemo-sim-kit48: agora CMB+tSZ [2025.02] take2 -- RA and DEC are unlensed +# nemo-sim-kit61: agora CMB+tSZ+kSZ [2026.03] +# nemo-sim-kit62: agora CMB+tSZ+CIB [2026.03] +# nemo-sim-kit63: agora CMB+tSZ+kSZ+CIB [2026.03] + +simsuite_version: v11 +# v11: scripts are tidied and merged [2026.06] + +ivar090: /data5/act/maps/dr5/act_planck_dr5.01_s08s18_AA_f090_night_ivar.fits +ivar150: /data5/act/maps/dr5/act_planck_dr5.01_s08s18_AA_f150_night_ivar.fits + + +agora_true_kappa: kappa/kappa_alm_recon_agora_phiNG_phi1_seed1.fits +agora_lensed_cmb: cmb/len/tqu1/agora_tqu1_phiNG_seed1_lmax16000_nside8192_interp1.6_method1_pol_1_lensedmap.fits +agora_tsz: tsz/unl/agora_utszNG_bahamas80_bnd_unb_1.0e+12_1.0e+18.fits +agora_ksz: ksz/unl/agora_ukszNG_bahamas80_bnd_unb_1.0e+12_1.0e+18.fits +agora_cib090: cib/uK/len/act/agora_len_mag_cibmap_act_90ghz_uk.fits +agora_cib150: cib/uK/len/act/agora_len_mag_cibmap_act_150ghz_uk.fits + + +# model_f090: nemo-sim-kit48/agora_cmb_tsz/clusterModelMaps/clusterModelMap_f090.fits +# model_f150: nemo-sim-kit48/agora_cmb_tsz/clusterModelMaps/clusterModelMap_f150.fits +# model_ksz_f090: nemo-sim-kit61/20260320_cmb_tsz_ksz/20260320_cmb_tsz_ksz/clusterModelMaps/clusterModelMap_f090.fits +# model_ksz_f150: nemo-sim-kit61/20260320_cmb_tsz_ksz/20260320_cmb_tsz_ksz/clusterModelMaps/clusterModelMap_f150.fits +# model_cib_f090: nemo-sim-kit62/20260320_cmb_tsz_cib/20260320_cmb_tsz_cib/clusterModelMaps/clusterModelMap_f090.fits +# model_cib_f150: nemo-sim-kit62/20260320_cmb_tsz_cib/20260320_cmb_tsz_cib/clusterModelMaps/clusterModelMap_f150.fits +# model_ksz_cib_f090: nemo-sim-kit63/20260320_cmb_tsz_ksz_cib/20260320_cmb_tsz_ksz_cib/clusterModelMaps/clusterModelMap_f090.fits +# model_ksz_cib_f150: nemo-sim-kit63/20260320_cmb_tsz_ksz_cib/20260320_cmb_tsz_ksz_cib/clusterModelMaps/clusterModelMap_f150.fits diff --git a/make_sim_suite.py b/make_sim_suite.py index 4957942..0c1ce18 100644 --- a/make_sim_suite.py +++ b/make_sim_suite.py @@ -1,7 +1,4 @@ import time as t -import matplotlib -matplotlib.use("Agg") -import matplotlib.pyplot as plt import numpy as np from pixell import enmap, reproject, utils, curvedsky, wcsutils, bunch from orphics import maps, io @@ -16,9 +13,7 @@ print("Paths: ", paths) - - -# RUN SETTING ------------------------------------------------------------------ +# RUN SETTING -------------------------------------------------------------------------------- parser = argparse.ArgumentParser() parser.add_argument( @@ -28,350 +23,315 @@ "which_sim", type=str, help="Choose the sim e.g. websky or sehgal or agora." ) parser.add_argument( - "--high-accuracy", action="store_true", help="If using a high accuracy websky map, choose this option." + "footprint", type=str, help="Choose the footprint e.g. full or act. The fullsky map is needed for halo sample." ) -args = parser.parse_args() +args = parser.parse_args() + output_path = paths.simsuite_path save_name = args.save_name -if args.which_sim == "websky": - print(" ::: producing maps for WEBSKY sim") - save_name = "websky_" + save_name - sim_path = paths.websky_sim_path -elif args.which_sim == "sehgal": - print(" ::: producing maps for SEHGAL sim") - save_name = "sehgal_" + save_name - sim_path = paths.sehgal_sim_path -elif args.which_sim == "agora": +if args.which_sim == "agora": print(" ::: producing maps for AGORA sim") save_name = "agora_" + save_name sim_path = paths.agora_sim_path save_dir = f"{output_path}/{save_name}/" io.mkdir(f"{save_dir}") -print(" ::: saving to", save_dir) - - +print(" ::: saving to", save_dir) - - -# SIM SETTING ------------------------------------------------------------------ - -freq_sz = 150 -print(" ::: map frequency at %d GHz" %freq_sz) +# SIM SETTING -------------------------------------------------------------------------------- px = 0.5 -fwhm = 1.5 -nlevel = 15.0 - +fwhm090 = 2.2 +fwhm150 = 1.5 +fwhm_plc = 5.0 +nlevel = 15.0 +nlevel_plc = 35.0 +# flux per pixel [mJy] cut for CIB map +fluxcut_val = 5.0 -# BEAM CONVOLUTION ------------------------------------------------------------- +# BEAM CONVOLUTION --------------------------------------------------------------------------- -def apply_beam(imap): +def apply_beam(imap, ifwhm): # map2alm of the maps, almxfl(alm, beam_1d) to convolve with beam, alm2map to convert back to map if args.which_sim == "websky" or args.which_sim == "agora": nside = 8192 elif args.which_sim == "sehgal": nside = 4096 alm_lmax = nside * 3 - bfunc = lambda x: maps.gauss_beam(fwhm, x) + bfunc = lambda x: maps.gauss_beam(ifwhm, x) imap_alm = curvedsky.map2alm(imap, lmax=alm_lmax) - beam_convoloved_alm = curvedsky.almxfl(imap_alm, bfunc) - return curvedsky.alm2map(beam_convoloved_alm, enmap.empty(imap.shape, imap.wcs)) - - - -# PREPARING MAPS --------------------------------------------------------------- - -# paths for reprojected maps -if args.which_sim == "websky": - if args.high_accuracy: - r_lmap = sim_path + paths.websky_dlensed_reproj - r_kmap = sim_path + paths.websky_kappa4p5_reproj - r_tszmap = sim_path + paths.websky_dtsz_reproj - r_kszmap = sim_path + paths.websky_dksz_reproj - r_cib093map = sim_path + paths.websky_dcib093_reproj - r_cib145map = sim_path + paths.websky_dcib145_reproj - else: - r_lmap = sim_path + paths.websky_cmb_reproj - r_kmap = sim_path + paths.websky_kappa_reproj - r_tszmap = sim_path + paths.websky_tsz_reproj -elif args.which_sim == "sehgal": - r_lmap = sim_path + paths.sehgal_cmb_reproj - r_kmap = sim_path + paths.sehgal_kappa_reproj - r_tszmap = sim_path + paths.sehgal_tsz_reproj -elif args.which_sim == "agora": - r_lmap = sim_path + paths.agora_cmb_reproj - r_kmap = sim_path + paths.agora_kappa_reproj - r_tszmap = sim_path + paths.agora_tsz_reproj - - -# reading lensed cmb map -try: - lmap = enmap.read_map(r_lmap, delayed=False) - print(" ::: reading reprojected lensed cmb map:", r_lmap) - -except: - shape, wcs = enmap.fullsky_geometry(res=px*utils.arcmin, proj="car") - - if args.which_sim == "websky": - if args.high_accuracy: - ifile = f"{paths.mat_path}/dlensed.fits" - lmap = enmap.read_map(ifile) - print(" ::: reading lensed cmb map:", ifile) - else: - ifile = f"{sim_path}lensed_alm.fits" - alm = np.complex128(hp.read_alm(ifile, hdu=(1, 2, 3))) - lmap = curvedsky.alm2map(alm[0,:], enmap.empty(shape, wcs, dtype=np.float64)) - print(" ::: reading lensed alm map and converting to lensed cmb map:", ifile) - - elif args.which_sim == "sehgal": - ifile = f"{sim_path}Sehgalsimparams_healpix_4096_KappaeffLSStoCMBfullsky_phi_SimLens_Tsynfastnopell_fast_lmax8000_nside4096_interp2.5_method1_1_lensed_map.fits" - hmap = hp.read_map(ifile).astype(np.float64) - lmap = reproject.healpix2map(hmap, shape, wcs) - print(" ::: reading lensed cmb map:", ifile) - - elif args.which_sim == "agora": - ifile = f"{sim_path}cmb/len/tqu1/agora_tqu1_phiNG_seed1_lmax16000_nside8192_interp1.6_method1_pol_1_lensedmap.fits" - hmap,_,_ = hp.read_map(ifile, field=[0,1,2]).astype(np.float64) - lmap = reproject.healpix2map(hmap, shape, wcs) - print(" ::: reading lensed cmb map:", ifile) - - enmap.write_map(r_lmap, lmap) -print("lmap", lmap.shape) -shape, wcs = lmap.shape, lmap.wcs - - -# reading true kappa map -try: - kmap = enmap.read_map(r_kmap, delayed=False) - print(" ::: reading reprojected true kappa map:", r_kmap) - -except: - if args.which_sim == "websky": - if args.high_accuracy: - ifile = f"{sim_path}kap_lt4.5.fits" # CMB lensing convergence from z<4.5 from halo+field websky - else: - ifile = f"{sim_path}kap.fits" # CMB lensing convergence from 0 uK_CMB - # i.e., you would multiply a map in 1e-26 Jy/sr by this factor to get an output map in uK_CMB - nu *= 1e9 - X = hplanck*nu/(kboltz*tcmb) - dBnudT = (2.*hplanck*nu**3.)/clight**2. * (np.exp(X))/(np.exp(X)-1.)**2. * X/tcmb_uK * 1e26 - return 1./dBnudT +# PREPARING TRUE KAPPA ----------------------------------------------------------------------- - # frequency for CIB - if freq_sz == 150: freq_cib = 145 - elif freq_sz == 90: freq_cib = 93 +if args.which_sim == "agora": + ifile = f"{sim_path}{paths.agora_true_kappa}" + hmap = hp.read_alm(ifile).astype(np.complex128) + kmap = curvedsky.alm2map(hmap, enmap.empty(shape, wcs, dtype=np.float64)) - cibmap = ItoDeltaT(freq_cib) * cibmap_i + print(" ::: reading and reprojecting true kappa map:", ifile) + print("kmap", kmap.shape) - print(" ::: converting cib map at %d GHz" %freq_cib) - enmap.write_map(r_cibmap, cibmap) -print("cibmap", cibmap.shape) -assert wcsutils.equal(kmap.wcs, lmap.wcs) -assert wcsutils.equal(kmap.wcs, tszmap.wcs) -assert wcsutils.equal(kmap.wcs, kszmap.wcs) -assert wcsutils.equal(kmap.wcs, cibmap.wcs) +# PREPARING LENSED CMB MAP ------------------------------------------------------------------- -scmb = lmap -scmb_cib = lmap + cibmap -scmb_ksz_cib = lmap + kszmap + cibmap +if args.which_sim == "agora": + ifile = f"{sim_path}{paths.agora_lensed_cmb}" + hmap,_,_ = hp.read_map(ifile, field=[0,1,2]).astype(np.float64) + lmap = reproject.healpix2map(hmap, shape, wcs) -scmb_tsz = lmap + tszmap -scmb_tsz_cib = lmap + tszmap + cibmap -scmb_tsz_ksz_cib = lmap + tszmap + kszmap + cibmap + print(" ::: reading and reprojecting lensed cmb map:", ifile) + print("lmap", lmap.shape) -print(" ::: SIGNAL maps are ready!") -white_noise = maps.white_noise(lmap.shape, lmap.wcs, noise_muK_arcmin=nlevel) +# PREPARING TSZ MAP (UNLENSED) --------------------------------------------------------------- -print(" ::: applying beam and adding white noise of %.1f uK" %nlevel) +if args.which_sim == "agora": + ifile = f"{sim_path}{paths.agora_tsz}" + hmap = hp.read_map(ifile).astype(np.float64) + ymap = reproject.healpix2map(hmap, shape, wcs) -bcmb = apply_beam(scmb) -bcmb_ksz_cib = apply_beam(scmb_ksz_cib) + print(" ::: reading and reprojecting ymap:", ifile) + print("ymap", ymap.shape) -ocmb = apply_beam(scmb) + white_noise -ocmb_cib = apply_beam(scmb_cib) + white_noise -ocmb_ksz_cib = apply_beam(scmb_ksz_cib) + white_noise +# convert compton-y to delta-T (in uK) +tcmb = 2.726 +tcmb_uK = tcmb * 1e6 #micro-Kelvin +H_cgs = 6.62608e-27 +K_cgs = 1.3806488e-16 -ocmb_tsz = apply_beam(scmb_tsz) + white_noise -ocmb_tsz_cib = apply_beam(scmb_tsz_cib) + white_noise -ocmb_tsz_ksz_cib = apply_beam(scmb_tsz_ksz_cib) + white_noise +def fnu(nu): + """ + nu in GHz + tcmb in Kelvin + """ + mu = H_cgs*(1e9*nu)/(K_cgs*tcmb) + ans = mu/np.tanh(old_div(mu,2.0)) - 4.0 + return ans -print(" ::: OBSERVED maps are ready!") +def y_to_tsz(y_map, freq): + # convert compton-y map to tSZ temperature map + print(f" ::: converting ymap to tsz map at {freq:.0f} GHz") + return fnu(freq) * y_map * tcmb_uK -save_scmb = f"{save_dir}scmb.fits" -save_scmb_cib = f"{save_dir}scmb_cib.fits" -save_scmb_ksz_cib = f"{save_dir}scmb_ksz_cib.fits" +tszmap090 = y_to_tsz(ymap, 90.0) +tszmap150 = y_to_tsz(ymap, 150.0) -save_scmb_tsz = f"{save_dir}scmb_tsz.fits" -save_scmb_tsz_cib = f"{save_dir}scmb_tsz_cib.fits" -save_scmb_tsz_ksz_cib = f"{save_dir}scmb_tsz_ksz_cib.fits" -save_ocmb = f"{save_dir}ocmb.fits" -save_ocmb_cib = f"{save_dir}ocmb_cib.fits" -save_ocmb_ksz_cib = f"{save_dir}ocmb_ksz_cib.fits" +# PREPARING KSZ MAP (UNLENSED) --------------------------------------------------------------- -save_ocmb_tsz = f"{save_dir}ocmb_tsz.fits" -save_ocmb_tsz_cib = f"{save_dir}ocmb_tsz_cib.fits" -save_ocmb_tsz_ksz_cib = f"{save_dir}ocmb_tsz_ksz_cib.fits" +if args.which_sim == "agora": + ifile = f"{sim_path}{paths.agora_ksz}" + hmap = hp.read_map(ifile).astype(np.float64) + kszmap = reproject.healpix2map(hmap, shape, wcs) -save_bcmb = f"{save_dir}bcmb.fits" -save_bcmb_ksz_cib = f"{save_dir}bcmb_ksz_cib.fits" + print(" ::: reading and reprojecting ksz map:", ifile) + print("kszmap", kszmap.shape) -print(" ::: saving all maps") -enmap.write_map(save_scmb, scmb) -enmap.write_map(save_scmb_cib, scmb_cib) -enmap.write_map(save_scmb_ksz_cib, scmb_ksz_cib) -enmap.write_map(save_scmb_tsz, scmb_tsz) -enmap.write_map(save_scmb_tsz_cib, scmb_tsz_cib) -enmap.write_map(save_scmb_tsz_ksz_cib, scmb_tsz_ksz_cib) +# PREPARING CIB MAP (LENSED; can't find the unlensed map T_T) --------------------------------- +# flux cut routine is copied from Karen's script for now! -enmap.write_map(save_ocmb, ocmb) -enmap.write_map(save_ocmb_cib, ocmb_cib) -enmap.write_map(save_ocmb_ksz_cib, ocmb_ksz_cib) +if args.which_sim == "agora": -enmap.write_map(save_ocmb_tsz, ocmb_tsz) -enmap.write_map(save_ocmb_tsz_cib, ocmb_tsz_cib) -enmap.write_map(save_ocmb_tsz_ksz_cib, ocmb_tsz_ksz_cib) + def flux_density_to_temp(freq_GHz): + # get factor for converting delta flux density in [MJy/sr] to delta T in CMB units [uK] + freq = float(freq_GHz) + x = freq / 56.8 + return (1.05e3 * (np.exp(x)-1)**2 * + np.exp(-x) * (freq / 100)**-4) -enmap.write_map(save_bcmb, bcmb) -enmap.write_map(save_bcmb_ksz_cib, bcmb_ksz_cib) + def fluxcut_to_cib_map(freq, fluxcut_mJy): + if freq == 90: + ifile = f"{sim_path}{paths.agora_cib090}" # in uK + elif freq == 150: + ifile = f"{sim_path}{paths.agora_cib150}" # in uK + hmap = hp.read_map(ifile).astype(np.float64) + print(f"\n{freq} GHz") + print("before cut:", hmap.min(), hmap.max(), hmap.mean()) + + nside = hp.get_nside(hmap) # extract the HEALPix resolution + pixel_solid_angle = hp.nside2pixarea(nside) # in steradian + lim = fluxcut_mJy * 1e-9 / pixel_solid_angle # [mJy] to [MJy/sr] + uK_lim = flux_density_to_temp(freq) * lim # [MJy/sr] to [uK] + + print(f"conversion factor: {flux_density_to_temp(freq):.2f}") + print(f"uK threshold: {uK_lim:.2f}") + print(f"number of pixels above threshold: {np.sum(hmap > uK_lim)} ({np.sum(hmap > uK_lim)/np.sum(hmap)*100:.6f} percent)") + print(f"implementing flux cut on map of {fluxcut_mJy} mJy per pixel") + + hmap[hmap > uK_lim]=0. + print("after cut:", hmap.min(), hmap.max(), hmap.mean()) + + return reproject.healpix2map(hmap, shape, wcs) + + cibmap090 = fluxcut_to_cib_map(90, fluxcut_val) + cibmap150 = fluxcut_to_cib_map(150, fluxcut_val) + + +for m in ( + kmap, + tszmap090, + tszmap150, + kszmap, + cibmap090, + cibmap150, +): + assert wcsutils.equal(lmap.wcs, m.wcs) + + +# ADDING MAPS -------------------------------------------------------------------------------- + +signal_maps = { + "cmb": lmap, + "cmb_ksz": lmap + kszmap, + "cmb_cib090": lmap + cibmap090, + "cmb_cib150": lmap + cibmap150, + "cmb_ksz_cib090": lmap + kszmap + cibmap090, + "cmb_ksz_cib150": lmap + kszmap + cibmap150, + "cmb_tsz090": lmap + tszmap090, + "cmb_tsz150": lmap + tszmap150, + "cmb_tsz090_ksz": lmap + tszmap090 + kszmap, + "cmb_tsz150_ksz": lmap + tszmap150 + kszmap, + "cmb_tsz090_cib": lmap + tszmap090 + cibmap090, + "cmb_tsz150_cib": lmap + tszmap150 + cibmap150, + "cmb_tsz090_ksz_cib": lmap + tszmap090 + kszmap + cibmap090, + "cmb_tsz150_ksz_cib": lmap + tszmap150 + kszmap + cibmap150, +} + +print(" ::: SIGNAL only maps are ready!") + + + +# APPLYING BEAM AND ADDING NOISE ------------------------------------------------------------- + +def make_observed_map_for_recon(smap, fwhm, noise_level, seed): + return ( + apply_beam(smap, fwhm) + + maps.white_noise( + shape, wcs, + noise_muK_arcmin=noise_level, + seed=seed, + ) + ) + +obs_map_details = { + "g_ocmb": ("cmb", fwhm_plc, nlevel_plc, 100), + "g_ocmb_ksz": ("cmb_ksz", fwhm_plc, nlevel_plc, 100), + "g_ocmb_cib": ("cmb_cib150", fwhm_plc, nlevel_plc, 100), + "g_ocmb_ksz_cib": ("cmb_ksz_cib150", fwhm_plc, nlevel_plc, 100), + "h_ocmb": ("cmb", fwhm150, nlevel, 101), + "h_ocmb_ksz": ("cmb_ksz", fwhm150, nlevel, 101), + "h_ocmb_cib090": ("cmb_cib090", fwhm090, nlevel, 102), + "h_ocmb_cib150": ("cmb_cib150", fwhm150, nlevel, 103), + "h_ocmb_ksz_cib090": ("cmb_ksz_cib090", fwhm090, nlevel, 102), + "h_ocmb_ksz_cib150": ("cmb_ksz_cib150", fwhm150, nlevel, 103), + "h_ocmb_tsz090": ("cmb_tsz090", fwhm090, nlevel, 102), + "h_ocmb_tsz150": ("cmb_tsz150", fwhm150, nlevel, 103), + "h_ocmb_tsz090_ksz": ("cmb_tsz090_ksz", fwhm090, nlevel, 102), + "h_ocmb_tsz150_ksz": ("cmb_tsz150_ksz", fwhm150, nlevel, 103), + "h_ocmb_tsz090_cib": ("cmb_tsz090_cib", fwhm090, nlevel, 102), + "h_ocmb_tsz150_cib": ("cmb_tsz150_cib", fwhm150, nlevel, 103), + "h_ocmb_tsz090_ksz_cib": ("cmb_tsz090_ksz_cib", fwhm090, nlevel, 102), + "h_ocmb_tsz150_ksz_cib": ("cmb_tsz150_ksz_cib", fwhm150, nlevel, 103), +} + +omaps = { + name: make_observed_map_for_recon( + signal_maps[signal_name], + fwhm, + noise_level, + seed, + ) + for name, (signal_name, fwhm, noise_level, seed) + in obs_map_details.items() +} + +print(" ::: applying corresponding beams and adding flatwhite noise") +print(" ::: OBSERVED maps for lensing reconstruction are ready!") + + +if args.footprint == "act": + + ivar090 = enmap.read_map(paths.ivar090) + ivar150 = enmap.read_map(paths.ivar150) + + def make_noise(ivar, seed): + noise = maps.modulated_noise_map(ivar, seed=seed)[0] + noise = np.nan_to_num(noise, posinf=0, neginf=0) + assert np.all(np.isfinite(noise)) + return noise + + def make_observed_map(sky_map, fwhm, noise_map): + return apply_beam(sky_map, fwhm) + noise_map + + freq_info = { + "090": (fwhm090, make_noise(ivar090, seed=1)), + "150": (fwhm150, make_noise(ivar150, seed=0)), + } + + dr5white_omap_names = [ + "cmb_tsz090", + "cmb_tsz150", + "cmb_tsz090_ksz", + "cmb_tsz150_ksz", + "cmb_tsz090_cib", + "cmb_tsz150_cib", + "cmb_tsz090_ksz_cib", + "cmb_tsz150_ksz_cib", + ] + + dr5white_omaps = { + f"h_ocmb{name[3:]}": make_observed_map( + signal_maps[name], + *freq_info["090" if "090" in name else "150"], + ) + for name in dr5white_omap_names + } + + print(" ::: applying corresponding beams and adding dr5white noise") + print(" ::: OBSERVED maps for model subtraction are ready!") + + + +# SAVING MAPS -------------------------------------------------------------------------------- + +suffix = "_fullsky" if args.footprint == "full" else "" +enmap.write_map(f"{save_dir}/true_kappa{suffix}.fits", kmap) + +for name, omap in omaps.items(): + enmap.write_map(f"{save_dir}/{name}{suffix}.fits", omap) + +if args.footprint == "act": + for name, dr5white_omap in dr5white_omaps.items(): + enmap.write_map(f"{save_dir}/{name}{suffix}_dr5white.fits", dr5white_omap) + +print(" ::: all maps are saved! yayyyyyy") elapsed = t.time() - start -print("\r ::: entire run took %.1f seconds" %elapsed) - - - - - - - - - - - - - - - - +print("\r ::: entire run took %.1f seconds" %elapsed) \ No newline at end of file diff --git a/recon_sim.py b/recon_sim.py index 1a704a2..bcd53a6 100644 --- a/recon_sim.py +++ b/recon_sim.py @@ -208,15 +208,15 @@ def bin(data, modrmap, bin_edges): ras, decs, zs, masses = np.loadtxt(cat, unpack=True) -elif args.which_cat == "tsz_true": - if args.full_sample: - cat = cat_path + f"{args.which_sim}_tsz_true.txt" - elif args.highsnr_sample: - cat = cat_path + f"{args.which_sim}_tsz_snr7_true.txt" - else: - cat = cat_path + f"{args.which_sim}_tsz_snr5p5_true.txt" # M200c, true mass, lensed coordinates for Agora +# elif args.which_cat == "tsz_true": +# if args.full_sample: +# cat = cat_path + f"{args.which_sim}_tsz_true.txt" +# elif args.highsnr_sample: +# cat = cat_path + f"{args.which_sim}_tsz_snr7_true.txt" +# else: +# cat = cat_path + f"{args.which_sim}_tsz_snr5p5_true.txt" # M200c, true mass, lensed coordinates for Agora - ras, decs, zs, masses = np.loadtxt(cat, unpack=True) +# ras, decs, zs, masses = np.loadtxt(cat, unpack=True) elif args.which_cat == "tsz": if args.which_sim == "websky": @@ -229,13 +229,15 @@ def bin(data, modrmap, bin_edges): zs = hdu[1].data["redshift"] else: if args.which_sim == "sehgal": cat = paths.sehgal_tsz_cat - elif args.which_sim == "agora": cat = paths.agora_tsz_cat + elif args.which_sim == "agora": cat = f"{paths.nsk_path}{paths.nemosim_version}/{paths.agora_tsz_cat}" + + print(cat) hdu = fits.open(cat) ras = hdu[1].data["RADeg"] decs = hdu[1].data["DECDeg"] - masses = hdu[1].data["M200c"] # 1e14 Msun #FIXME: make it as an option - # masses = hdu[1].data["M200m"] # 1e14 Msun ##### + masses = hdu[1].data["M200c"] # 1e14 Msun + # masses = hdu[1].data["M200m"] # 1e14 Msun snr = hdu[1].data["SNR"] # fixed_SNR zs = hdu[1].data["redshift"] @@ -263,7 +265,7 @@ def bin(data, modrmap, bin_edges): Nx = 100 * len(ras) - if args.which_cat == "tsz_true": args.which_cat = "tsz" + # if args.which_cat == "tsz_true": args.which_cat = "tsz" # load random catalogue - created by mapcat.py + randcat.py cat = cat_path + f"{args.which_sim}_{args.which_cat}_randoms.txt" ras, decs = np.loadtxt(cat, unpack=True) @@ -280,19 +282,8 @@ def bin(data, modrmap, bin_edges): # READING MAPS ----------------------------------------------------------------- -if args.which_sim == "websky": - if args.high_accuracy: - true = paths.websky_sim_path + paths.websky_kappa4p5_reproj - else: - true = paths.websky_sim_path + paths.websky_kappa_reproj -elif args.which_sim == "sehgal": - if args.high_accuracy: - true = paths.sehgal_sim_path + paths.sehgal_dkappa_reproj - else: - true = paths.sehgal_sim_path + paths.sehgal_kappa_reproj -elif args.which_sim == "agora": - true = paths.agora_sim_path + paths.agora_kappa_reproj - +if args.which_sim == "agora": + true = f"{simsuite_path}true_kappa.fits" print(" ::: preparing for OBSERVED maps") @@ -308,21 +299,23 @@ def bin(data, modrmap, bin_edges): h_cmb_map = enmap.read_map(h_cmb, delayed=False) g_cmb_map = enmap.read_map(g_cmb, delayed=False) + if args.cmb_tsz: - cmb_tsz150 = f"{simsuite_path}h_ocmb_tsz150.fits" cmb_tsz090 = f"{simsuite_path}h_ocmb_tsz090.fits" + cmb_tsz150 = f"{simsuite_path}h_ocmb_tsz150.fits" - print(" ::: reading lensed cmb + tsz map at 150GHz:", cmb_tsz150) print(" ::: reading lensed cmb + tsz map at 90GHz:", cmb_tsz090) + print(" ::: reading lensed cmb + tsz map at 150GHz:", cmb_tsz150) - cmb_tsz_map150 = enmap.read_map(cmb_tsz150, delayed=False) cmb_tsz_map090 = enmap.read_map(cmb_tsz090, delayed=False) + cmb_tsz_map150 = enmap.read_map(cmb_tsz150, delayed=False) -if args.cmb_msub: - cmb_msub090 = f"{simsuite_path}h_ocmb_msub090_test.fits" - cmb_msub150 = f"{simsuite_path}h_ocmb_msub150_test.fits" +if args.cmb_msub: + + cmb_msub090 = f"{simsuite_path}h_ocmb_tsz090_msub.fits" + cmb_msub150 = f"{simsuite_path}h_ocmb_tsz150_msub.fits" print(" ::: reading the model subtracted maps!", cmb_msub090) print(" ::: reading the model subtracted maps!", cmb_msub150) @@ -330,6 +323,7 @@ def bin(data, modrmap, bin_edges): cmb_msub_map090 = enmap.read_map(cmb_msub090, delayed=False) cmb_msub_map150 = enmap.read_map(cmb_msub150, delayed=False) + if args.cmb_ksz: h_cmb_ksz = f"{simsuite_path}h_ocmb_ksz.fits" @@ -341,11 +335,12 @@ def bin(data, modrmap, bin_edges): h_cmb_ksz_map = enmap.read_map(h_cmb_ksz, delayed=False) g_cmb_ksz_map = enmap.read_map(g_cmb_ksz, delayed=False) + if args.cmb_cib: - h_cmb_cib090 = f"{simsuite_path}h_ocmb_cib090_5.fits" - h_cmb_cib150 = f"{simsuite_path}h_ocmb_cib150_5.fits" - g_cmb_cib = f"{simsuite_path}g_ocmb_cib5.fits" + h_cmb_cib090 = f"{simsuite_path}h_ocmb_cib090.fits" + h_cmb_cib150 = f"{simsuite_path}h_ocmb_cib150.fits" + g_cmb_cib = f"{simsuite_path}g_ocmb_cib.fits" print(" ::: reading lensed cmb + cib090 map for hres:", h_cmb_cib090) print(" ::: reading lensed cmb + cib150 map for hres:", h_cmb_cib150) @@ -370,11 +365,12 @@ def bin(data, modrmap, bin_edges): h_cmb_ksz_tsz150_map = enmap.read_map(h_cmb_ksz_tsz150, delayed=False) g_cmb_ksz_map = enmap.read_map(g_cmb_ksz, delayed=False) + if args.cmb_cib_tsz: - h_cmb_cib_tsz090 = f"{simsuite_path}h_ocmb_tsz090_cib5.fits" - h_cmb_cib_tsz150 = f"{simsuite_path}h_ocmb_tsz150_cib5.fits" - g_cmb_cib = f"{simsuite_path}g_ocmb_cib5.fits" + h_cmb_cib_tsz090 = f"{simsuite_path}h_ocmb_tsz090_cib.fits" + h_cmb_cib_tsz150 = f"{simsuite_path}h_ocmb_tsz150_cib.fits" + g_cmb_cib = f"{simsuite_path}g_ocmb_cib.fits" print(" ::: reading lensed cmb + cib + tsz090 map for hres:", h_cmb_cib_tsz090) print(" ::: reading lensed cmb + cib + tsz150 map for hres:", h_cmb_cib_tsz150) @@ -384,14 +380,14 @@ def bin(data, modrmap, bin_edges): h_cmb_cib_tsz150_map = enmap.read_map(h_cmb_cib_tsz150, delayed=False) g_cmb_cib_map = enmap.read_map(g_cmb_cib, delayed=False) + if args.cmb_ksz_msub: - # cmb_msub090 = f"{simsuite_path}h_ocmb_msub090_ksz_updated.fits" - # cmb_msub150 = f"{simsuite_path}h_ocmb_msub150_ksz_updated.fits" - cmb_msub090 = f"{simsuite_path}/../../maps/agora_h_ocmb_tsz_ksz_f090_msub.fits" # new model subtraction method - cmb_msub150 = f"{simsuite_path}/../../maps/agora_h_ocmb_tsz_ksz_f150_msub.fits" + # initially used the model map from nemo + # currently using a new model map generated from tsz_subtract.py + cmb_msub090 = f"{simsuite_path}h_ocmb_tsz090_ksz_msub.fits" + cmb_msub150 = f"{simsuite_path}h_ocmb_tsz150_ksz_msub.fits" g_cmb_ksz = f"{simsuite_path}g_ocmb_ksz.fits" - print(" ::: reading the model subtracted maps!", cmb_msub090) print(" ::: reading the model subtracted maps!", cmb_msub150) @@ -402,12 +398,11 @@ def bin(data, modrmap, bin_edges): g_cmb_ksz_map = enmap.read_map(g_cmb_ksz, delayed=False) - if args.cmb_cib_msub: - cmb_msub090 = f"{simsuite_path}h_ocmb_msub090_cib5_updated.fits" - cmb_msub150 = f"{simsuite_path}h_ocmb_msub150_cib5_updated.fits" - g_cmb_cib = f"{simsuite_path}g_ocmb_cib5.fits" + cmb_msub090 = f"{paths.msub_cib_f090}" + cmb_msub150 = f"{paths.msub_cib_f150}" + g_cmb_cib = f"{simsuite_path}g_ocmb_cib.fits" print(" ::: reading the model subtracted maps!", cmb_msub090) print(" ::: reading the model subtracted maps!", cmb_msub150) @@ -418,13 +413,11 @@ def bin(data, modrmap, bin_edges): g_cmb_cib_map = enmap.read_map(g_cmb_cib, delayed=False) - - if args.cmb_ksz_cib: - h_cmb_ksz_cib090 = f"{simsuite_path}h_ocmb_ksz_cib090_5.fits" - h_cmb_ksz_cib150 = f"{simsuite_path}h_ocmb_ksz_cib150_5.fits" - g_cmb_ksz_cib = f"{simsuite_path}g_ocmb_ksz_cib5.fits" + h_cmb_ksz_cib090 = f"{simsuite_path}h_ocmb_ksz_cib090.fits" + h_cmb_ksz_cib150 = f"{simsuite_path}h_ocmb_ksz_cib150.fits" + g_cmb_ksz_cib = f"{simsuite_path}g_ocmb_ksz_cib.fits" print(" ::: reading lensed cmb + ksz + cib090 map for hres:", h_cmb_ksz_cib090) print(" ::: reading lensed cmb + ksz + cib150 map for hres:", h_cmb_ksz_cib150) @@ -434,11 +427,12 @@ def bin(data, modrmap, bin_edges): h_cmb_ksz_cib150_map = enmap.read_map(h_cmb_ksz_cib150, delayed=False) g_cmb_ksz_cib_map = enmap.read_map(g_cmb_ksz_cib, delayed=False) + if args.cmb_ksz_cib_tsz: - h_cmb_tsz_ksz_cib090 = f"{simsuite_path}h_ocmb_tsz090_ksz_cib5.fits" - h_cmb_tsz_ksz_cib150 = f"{simsuite_path}h_ocmb_tsz150_ksz_cib5.fits" - g_cmb_ksz_cib = f"{simsuite_path}g_ocmb_ksz_cib5.fits" + h_cmb_tsz_ksz_cib090 = f"{simsuite_path}h_ocmb_tsz090_ksz_cib.fits" + h_cmb_tsz_ksz_cib150 = f"{simsuite_path}h_ocmb_tsz150_ksz_cib.fits" + g_cmb_ksz_cib = f"{simsuite_path}g_ocmb_ksz_cib.fits" print(" ::: reading lensed cmb + tsz090 + ksz + cib090 map for hres:", h_cmb_tsz_ksz_cib090) print(" ::: reading lensed cmb + tsz150 + ksz + cib150 map for hres:", h_cmb_tsz_ksz_cib150) @@ -448,12 +442,12 @@ def bin(data, modrmap, bin_edges): h_cmb_tsz_ksz_cib150_map = enmap.read_map(h_cmb_tsz_ksz_cib150, delayed=False) g_cmb_ksz_cib_map = enmap.read_map(g_cmb_ksz_cib, delayed=False) -if args.cmb_ksz_cib_msub: - cmb_msub090 = f"{simsuite_path}h_ocmb_msub090_ksz_cib5_updated.fits" - cmb_msub150 = f"{simsuite_path}h_ocmb_msub150_ksz_cib5_updated.fits" - g_cmb_ksz_cib = f"{simsuite_path}g_ocmb_ksz_cib5.fits" +if args.cmb_ksz_cib_msub: + cmb_msub090 = f"{paths.msub_ksz_cib_f090}" + cmb_msub150 = f"{paths.msub_ksz_cib_f150}" + g_cmb_ksz_cib = f"{simsuite_path}g_ocmb_ksz_cib.fits" print(" ::: reading the model subtracted maps!", cmb_msub090) print(" ::: reading the model subtracted maps!", cmb_msub150) @@ -464,7 +458,6 @@ def bin(data, modrmap, bin_edges): g_cmb_ksz_cib_map = enmap.read_map(g_cmb_ksz_cib, delayed=False) - print(" ::: maps are ready!") @@ -598,7 +591,7 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): res=px * utils.arcmin, proj="tan", oversample=2, - pixwin=True if args.high_accuracy else False # only True for maps made natively in CAR + pixwin=True if args.high_accuracy else False ) # lensed cmb for hres g_cmb = reproject.thumbnails( @@ -608,50 +601,52 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): res=px * utils.arcmin, proj="tan", oversample=2, - pixwin=True if args.high_accuracy else False # only True for maps made natively in CAR + pixwin=True if args.high_accuracy else False ) # lensed cmb for grad if args.cmb_tsz: - h_cmb_tsz150 = reproject.thumbnails( - cmb_tsz_map150, + h_cmb_tsz090 = reproject.thumbnails( + cmb_tsz_map090, coords, r=maxr, res=px * utils.arcmin, proj="tan", oversample=2, pixwin=False - ) # lensed cmb + tsz150 + ) # lensed cmb + tsz090 - h_cmb_tsz090 = reproject.thumbnails( - cmb_tsz_map090, + h_cmb_tsz150 = reproject.thumbnails( + cmb_tsz_map150, coords, r=maxr, res=px * utils.arcmin, proj="tan", oversample=2, pixwin=False - ) # lensed cmb + tsz090 + ) # lensed cmb + tsz150 + if args.cmb_msub: - h_cmb_msub150 = reproject.thumbnails( - cmb_msub_map150, + h_cmb_msub090 = reproject.thumbnails( + cmb_msub_map090, coords, r=maxr, res=px * utils.arcmin, proj="tan", oversample=2, pixwin=False - ) # lensed cmb + tsz150 - model150 + ) # lensed cmb + tsz090 - model090 - h_cmb_msub090 = reproject.thumbnails( - cmb_msub_map090, + h_cmb_msub150 = reproject.thumbnails( + cmb_msub_map150, coords, r=maxr, res=px * utils.arcmin, proj="tan", oversample=2, pixwin=False - ) # lensed cmb + tsz090 - model090 + ) # lensed cmb + tsz150 - model150 + if args.cmb_ksz: h_cmb_ksz = reproject.thumbnails( @@ -705,9 +700,6 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): pixwin=False ) # lensed cmb + cib150 for grad - - - if args.cmb_ksz_tsz: h_cmb_ksz_tsz090 = reproject.thumbnails( h_cmb_ksz_tsz090_map, @@ -748,7 +740,7 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): proj="tan", oversample=2, pixwin=False - ) # lensed cmb + cib + tsz090 for hres + ) # lensed cmb + cib090 + tsz090 for hres h_cmb_cib_tsz150 = reproject.thumbnails( h_cmb_cib_tsz150_map, @@ -758,7 +750,7 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): proj="tan", oversample=2, pixwin=False - ) # lensed cmb + cib + tsz150 for hres + ) # lensed cmb + cib150 + tsz150 for hres g_cmb_cib = reproject.thumbnails( g_cmb_cib_map, @@ -768,7 +760,7 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): proj="tan", oversample=2, pixwin=False - ) # lensed cmb + cib for grad + ) # lensed cmb + cib150 for grad if args.cmb_ksz_msub: h_cmb_ksz_msub090 = reproject.thumbnails( @@ -832,11 +824,6 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): pixwin=False ) # lensed cmb + cib for grad - - - - - if args.cmb_ksz_cib: h_cmb_ksz_cib090 = reproject.thumbnails( h_cmb_ksz_cib090_map, @@ -979,7 +966,6 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): assert wcsutils.equal(kstamp.wcs, h_cmb_cib_msub090.wcs) assert wcsutils.equal(kstamp.wcs, g_cmb_cib.wcs) - elif args.cmb_ksz_cib: assert wcsutils.equal(kstamp.wcs, h_cmb_ksz_cib150.wcs) assert wcsutils.equal(kstamp.wcs, h_cmb_ksz_cib090.wcs) @@ -1082,7 +1068,6 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): hres150 = h_cmb_cib150 hres90 = h_cmb_cib090 g_cmb = g_cmb_cib - elif args.cmb_ksz_tsz: hres150 = h_cmb_ksz_tsz150 hres90 = h_cmb_ksz_tsz090 @@ -1099,8 +1084,6 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): hres150 = h_cmb_cib_msub150 hres90 = h_cmb_cib_msub090 g_cmb = g_cmb_cib - - elif args.cmb_ksz_cib: hres150 = h_cmb_ksz_cib150 hres90 = h_cmb_ksz_cib090 @@ -1133,9 +1116,6 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): tapered_hres = hres * taper - - - s.add_to_stack("grad2d_before", g_cmb) g_filtered = maps.filter_map(g_cmb, ymask) s.add_to_stack("grad2d_before_filtered", g_filtered) @@ -1158,7 +1138,6 @@ def ilc(modlmap, m1, m2, p11, p22, p12, b1, b2): # taper stamp tapered_grad = g_cmb * taper - s.add_to_stack("grad2d_after", g_cmb) g_filtered = maps.filter_map(g_cmb, ymask) s.add_to_stack("grad2d_after_filtered", g_filtered)