From 0110261cee172e3c3c861018e055aa11a449f455 Mon Sep 17 00:00:00 2001 From: Neha Joshi Date: Thu, 8 May 2025 15:41:15 -0400 Subject: [PATCH 01/10] fetched from main branch --- kapstack.py | 54 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 54 insertions(+) create mode 100644 kapstack.py diff --git a/kapstack.py b/kapstack.py new file mode 100644 index 0000000..ba285c0 --- /dev/null +++ b/kapstack.py @@ -0,0 +1,54 @@ +import os,sys +import numpy as np +import healpy as hp +from orphics import maps,io,interfaces as ints,lensing,stats +from pixell import utils,bench,enmap + +rstamp = 30.0 * utils.arcmin +res = 0.25 * utils.arcmin + +with bench.show("load kappa"): + kappa = hp.read_map("/data5/sims/agora_sims/cmbkappa/raytrace16384_ip20_cmbkappa.zs1.kg1g2_highzadded_lowLcorrected.fits") + +z_min = 0.5 +z_max = 0.54 +mmin = 3e14 +mmax = 3.4e14 + +oras,odecs,ozs,oms = ints.get_agora_halos(z_min = z_min, z_max = z_max, + mass_min = mmin, mass_max = mmax, + massdef='m500') # Get agora halos + +print(f"Starting stack on {len(oras)}") +out = 0. +for i,(ora,odec) in enumerate(zip(oras,odecs)): + othumb = maps.thumbnail_healpix(kappa,np.asarray((odec,ora))*utils.degree,r=rstamp,res=res) # reproject directly from healpix + if i%100==0: print(f"Done {i+1} / {len(oras)}") + out = out + othumb + +out = out / len(oras) + +lmin = 200 +lmax = 3500 +rmin = 0. +rmax = 10*utils.arcmin +rwidth = 1.0 * utils.arcmin + +cents,kap_agora = lensing.filter_bin_kappa2d(out,lmin=lmin,lmax=lmax,rmin=0.,rmax=rmax,rwidth=rwidth) + + +z = ozs.mean() +ez = ozs.std() +m500c = oms.mean() +em500c = oms.std() + +thetas,kappa_1h,kappa_2h,tot_kappa,cents_t,b1d1h,b1d2h,b1d_t = lensing.kappa_nfw_profiley(mass=m500c,conc=None, + z=z,z_s=1100.,background='critical',delta=500,apply_filter=True, + lmin=lmin,lmax=lmax,res=res, + rstamp=rstamp,rmin=rmin,rmax=rmax,rwidth=rwidth) +pl = io.Plotter(xyscale='linlin',xlabel='$\\theta$ (arcmin)',ylabel='$\\kappa$') +pl.add(cents/utils.arcmin,kap_agora,ls='-',color='r',label=f'Agora stack; M={m500c/1e14:.1f} $\\pm$ {em500c/1e14:.1f} ; z={z:.1f} $\\pm$ {ez:.1f} ; ') +pl.add(cents_t/utils.arcmin,b1d_t,ls='--',color='k',label='profiley theory') +pl.hline(y=0) +pl.done('profile.png') +enmap.write_map("stack.fits",out) From ed1b66ec300a33fe008842a8e7e4b3103840766a Mon Sep 17 00:00:00 2001 From: Neha Joshi Date: Fri, 15 Aug 2025 16:12:38 -0400 Subject: [PATCH 02/10] updated flatsky sim case to make user's choice of template for lensing recon --- analysis.py | 27 ++++++++++++++++++++------- 1 file changed, 20 insertions(+), 7 deletions(-) diff --git a/analysis.py b/analysis.py index ea20e22..458d1da 100644 --- a/analysis.py +++ b/analysis.py @@ -39,11 +39,11 @@ def __init__(self, self.debug = debug - # Load template - 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_correlations + thetas_cross,kappa_1h_cross,kappa_2h_cross = np.loadtxt(c.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 @@ -53,7 +53,20 @@ def __init__(self, # will be slightly different from that constructed # directly - self.csim = olensing.FixedLens(thetas, kappa, width_deg=c.width_deg, pad_fact=1 if c.periodic else c.pad_fact) + # Make template profile for lensing + mass = c.mass + z=c.z + delta = c.delta + + thetas,kappa_1h,kappa_2h,_,_,_,_,_ = olensing.kappa_nfw_profiley(mass=mass,conc=None, + z=z,z_s=1100.,background='critical',delta=delta, R_off_Mpc = None, + apply_filter=False) + + kappa = kappa_1h + kappa_2h + self.thetas = thetas + self.template_kappa_1d = kappa + + self.csim = olensing.FixedLens(self.thetas, self.kappa, width_deg=c.width_deg, pad_fact=1 if c.periodic else c.pad_fact) _,_,dummy = self.csim.generate_sim(0) # FIXME: is this necessary self.shape, self.wcs = dummy.shape, dummy.wcs else: @@ -66,7 +79,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 = enmap.enmap(maps.interp(self.thetas_cross,self.kappa_cross)(self.modrmap),self.wcs) self.reconstructor = Recon(self.shape,self.wcs, From be0c554e56d520465e83cdff5b4fb36436267afb Mon Sep 17 00:00:00 2001 From: Neha Joshi Date: Wed, 20 Aug 2025 12:21:36 -0400 Subject: [PATCH 03/10] updating analysis.finish for separate reconstruction template --- analysis.py | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/analysis.py b/analysis.py index 458d1da..635e3db 100644 --- a/analysis.py +++ b/analysis.py @@ -43,7 +43,7 @@ def __init__(self, thetas_cross,kappa_1h_cross,kappa_2h_cross = np.loadtxt(c.template_file,unpack=True) kappa_cross = kappa_1h_cross + kappa_2h_cross self.thetas_cross = thetas_cross - self.template_kappa_1d_cross = kappa_cross + self.kappa_cross = kappa_cross if atype=='flatsky_sim': # If we are doing a padded high-res flat sim @@ -64,7 +64,7 @@ def __init__(self, kappa = kappa_1h + kappa_2h self.thetas = thetas - self.template_kappa_1d = kappa + self.kappa = kappa self.csim = olensing.FixedLens(self.thetas, self.kappa, width_deg=c.width_deg, pad_fact=1 if c.periodic else c.pad_fact) _,_,dummy = self.csim.generate_sim(0) # FIXME: is this necessary @@ -163,6 +163,7 @@ def get_recon(self,imapx,imapy=None,do_real=False,debug=False): # Cross-correlate with template cents,cp1d,cp2d = self.power(self.ktemplate,kmap) + self.template_cross_p1d = cp1d if not(np.all(np.isfinite(cp1d))): raise ValueError if debug and self.rank==0: io.plot_img((np.fft.fftshift(self.modlmap**2 * cp2d)),self._out('recon_template_cross_p2d.png')) @@ -251,7 +252,7 @@ def finish(self, save_transfer=False): lcents = self.lcents Lmin = self.Lmin Lmax = self.Lmax - p1d = self.template_auto_p1d + p1d = self.template_cross_p1d cmean = self.s.stats['cp1d']['mean'] cerr = self.s.stats['cp1d']['errmean'] ccov = self.s.stats['cp1d']['covmean'] @@ -317,13 +318,14 @@ def finish(self, save_transfer=False): pl._ax.set_ylim(-3e-5,5.1e-5) pl.done(self._out('recon_template_cross_p1d.png')) - pl = io.Plotter() - pl.add_err(rcents[rsel]/u.arcmin,rmean,yerr=rerr,marker='o') - pl.add(self.thetas/u.arcmin, self.template_kappa_1d) - pl._ax.set_xlim(0,10) - pl._ax.set_ylim(-0.01,rmean.max()*1.2) - pl.hline(y=0) - pl.done(self._out('recon_profile.png')) + if self.atype == 'flatsky_sim': + pl = io.Plotter() + pl.add_err(rcents[rsel]/u.arcmin,rmean,yerr=rerr,marker='o') + pl.add(self.thetas/u.arcmin, self.kappa) + pl._ax.set_xlim(0,10) + pl._ax.set_ylim(-0.01,rmean.max()*1.2) + pl.hline(y=0) + pl.done(self._out('recon_profile.png')) class Recon(object): From 19fda6feb3aa4042dae2e085ce096d10acc467e7 Mon Sep 17 00:00:00 2001 From: Neha Joshi Date: Wed, 20 Aug 2025 23:45:18 -0400 Subject: [PATCH 04/10] additional signal calc changes for template --- analysis.py | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/analysis.py b/analysis.py index 635e3db..5b584b7 100644 --- a/analysis.py +++ b/analysis.py @@ -69,6 +69,7 @@ def __init__(self, self.csim = olensing.FixedLens(self.thetas, self.kappa, width_deg=c.width_deg, pad_fact=1 if c.periodic else c.pad_fact) _,_,dummy = self.csim.generate_sim(0) # FIXME: is this necessary self.shape, self.wcs = dummy.shape, dummy.wcs + else: if c.periodic: raise ValueError self.thumbnail_r = c.width_deg * u.deg / 2. @@ -163,7 +164,6 @@ def get_recon(self,imapx,imapy=None,do_real=False,debug=False): # Cross-correlate with template cents,cp1d,cp2d = self.power(self.ktemplate,kmap) - self.template_cross_p1d = cp1d if not(np.all(np.isfinite(cp1d))): raise ValueError if debug and self.rank==0: io.plot_img((np.fft.fftshift(self.modlmap**2 * cp2d)),self._out('recon_template_cross_p2d.png')) @@ -252,7 +252,7 @@ def finish(self, save_transfer=False): lcents = self.lcents Lmin = self.Lmin Lmax = self.Lmax - p1d = self.template_cross_p1d + p1d = self.template_auto_p1d cmean = self.s.stats['cp1d']['mean'] cerr = self.s.stats['cp1d']['errmean'] ccov = self.s.stats['cp1d']['covmean'] @@ -274,7 +274,7 @@ def finish(self, save_transfer=False): xlabel='$\\Theta$ (arcmin)',ylabel='$\\Theta$ (arcmin)') sel = np.logical_and(lcents>Lmin,lcents Date: Wed, 20 Aug 2025 23:59:21 -0400 Subject: [PATCH 05/10] proper formatting --- choices.yaml | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/choices.yaml b/choices.yaml index b1efdb0..61e0b30 100644 --- a/choices.yaml +++ b/choices.yaml @@ -42,7 +42,10 @@ ybeam_noise_uk_arcmin: 10.0 # ================================================ # How the template and reconstruction are treated # ================================================ - +# This is the profile used for lensing reconstruction in the flatsky sim case +mass: 2.0e+14 +delta: 500 #'200' or '500' for critical mass (profiley will calculate m200c or m500c) +z: 1.0 theta_max_arc: 20.0 # All templates and reconstructions are masked beyond this radius so that we are not sensitive to the modeling there roll_arc: 20.0 # with a cosine taper that rolls of to zero in this distance # The above two decide what ctaper looks like (saved to ctaper.png) From bf954122fba6ced7b86222462de94eb5954cce12 Mon Sep 17 00:00:00 2001 From: Neha Joshi Date: Thu, 21 Aug 2025 12:32:14 -0400 Subject: [PATCH 06/10] fixed cross-corr for snr calculation --- analysis.py | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/analysis.py b/analysis.py index 5b584b7..f8e4f83 100644 --- a/analysis.py +++ b/analysis.py @@ -111,9 +111,11 @@ def __init__(self, self.Lmin = c.Lmin self.Lmax = c.Lmax + # Auto-correlation for template self.ktemplate = enmap.fft(self.template * self.ctaper,normalize='phys') cents,p1d,p2d = self.power(self.ktemplate,self.ktemplate) self.template_auto_p1d = p1d + if debug and self.rank==0: io.plot_img(self.ctaper,self._out('ctaper.png'),arc_width=c.width_deg*60., xlabel='$\\Theta$ (arcmin)',ylabel='$\\Theta$ (arcmin)') @@ -124,7 +126,12 @@ def __init__(self, pl.hline(y=0) pl.done(self._out('template_auto_p1d.png')) - + # Cross-correlation of template and lensing template + self.lensing_template = enmap.enmap(maps.interp(self.thetas,self.kappa)(self.modrmap),self.wcs) + self.klensing = enmap.fft(self.lensing_template * self.ctaper,normalize='phys') + cents, p1d_cross, p2d_cross = self.power(self.ktemplate, self.klensing) + self.template_cross_p1d = p1d_cross + self.fit_Kmask = maps.mask_kspace(self.shape, self.wcs, lmin=c.Lmin, lmax=c.Lmax) if self.rank==0: @@ -252,7 +259,9 @@ def finish(self, save_transfer=False): lcents = self.lcents Lmin = self.Lmin Lmax = self.Lmax - p1d = self.template_auto_p1d + if self.atype=='flatsky_sim': + p1d = self.template_cross_p1d + else: p1d = self.template_auto_p1d cmean = self.s.stats['cp1d']['mean'] cerr = self.s.stats['cp1d']['errmean'] ccov = self.s.stats['cp1d']['covmean'] @@ -274,7 +283,7 @@ def finish(self, save_transfer=False): xlabel='$\\Theta$ (arcmin)',ylabel='$\\Theta$ (arcmin)') sel = np.logical_and(lcents>Lmin,lcents Date: Thu, 21 Aug 2025 13:39:50 -0400 Subject: [PATCH 07/10] added option for unlensed flatsky sims --- analysis.py | 3 ++- choices.yaml | 3 +++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/analysis.py b/analysis.py index f8e4f83..4661026 100644 --- a/analysis.py +++ b/analysis.py @@ -66,7 +66,8 @@ def __init__(self, self.thetas = thetas self.kappa = kappa - self.csim = olensing.FixedLens(self.thetas, self.kappa, width_deg=c.width_deg, pad_fact=1 if c.periodic else c.pad_fact) + self.csim = olensing.FixedLens(self.thetas, self.kappa, width_deg=c.width_deg, pad_fact=1 if c.periodic else c.pad_fact, + skip_lensing=c.skip_lensing) _,_,dummy = self.csim.generate_sim(0) # FIXME: is this necessary self.shape, self.wcs = dummy.shape, dummy.wcs diff --git a/choices.yaml b/choices.yaml index 61e0b30..47dc8bd 100644 --- a/choices.yaml +++ b/choices.yaml @@ -42,6 +42,9 @@ ybeam_noise_uk_arcmin: 10.0 # ================================================ # How the template and reconstruction are treated # ================================================ +# Set whether lensing will be used in the flatsky sim case +skip_lensing: False + # This is the profile used for lensing reconstruction in the flatsky sim case mass: 2.0e+14 delta: 500 #'200' or '500' for critical mass (profiley will calculate m200c or m500c) From 1581e72b5976a6a10bdf2e37e9b3ff743532645e Mon Sep 17 00:00:00 2001 From: Neha Joshi Date: Thu, 21 Aug 2025 14:48:37 -0400 Subject: [PATCH 08/10] plot template cross and auto correlations --- analysis.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/analysis.py b/analysis.py index 4661026..1418149 100644 --- a/analysis.py +++ b/analysis.py @@ -133,6 +133,18 @@ def __init__(self, cents, p1d_cross, p2d_cross = self.power(self.ktemplate, self.klensing) self.template_cross_p1d = p1d_cross + if debug and self.rank==0: + cents, p1d_auto, p2d_auto = self.power(self.klensing, self.klensing) + + io.plot_img((np.fft.fftshift(self.modlmap**2 * p2d_cross)),self._out('template_cross_p2d.png')) + pl = io.Plotter(xyscale='linlin',xlabel='L',ylabel='L^2 C_L') + pl.add(cents,cents**2 * p1d_cross,marker='o',label='cross') + pl.add(cents,cents**2 * p1d_auto,marker='o',label='lensing template auto') + pl.add(cents,cents**2 * p1d,marker='o',label='correlation template auto') + pl.vline(x=500) + pl.hline(y=0) + pl.done(self._out('template_cross_p1d.png')) + self.fit_Kmask = maps.mask_kspace(self.shape, self.wcs, lmin=c.Lmin, lmax=c.Lmax) if self.rank==0: From a29fd7c336042eae87e7d9f81d79c9a4ce9eb12a Mon Sep 17 00:00:00 2001 From: Neha Joshi Date: Thu, 21 Aug 2025 15:20:53 -0400 Subject: [PATCH 09/10] kapstack.py --- kapstack.py | 96 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 96 insertions(+) diff --git a/kapstack.py b/kapstack.py index e69de29..28b7578 100644 --- a/kapstack.py +++ b/kapstack.py @@ -0,0 +1,96 @@ +import os,sys +import numpy as np +import healpy as hp +from orphics import maps,io,interfaces as ints,lensing,stats +from pixell import utils,bench,enmap,reproject + +rstamp = 30.1 * utils.arcmin +res = 0.1 * utils.arcmin +Nmax = None +#Nmax = 3 + +with bench.show("load kappa"): + kappa = hp.read_map("/data5/sims/agora_sims/cmbkappa/raytrace16384_ip20_cmbkappa.zs1.kg1g2_highzadded_lowLcorrected.fits") + npix = kappa.size + nside = hp.npix2nside(npix) + print("nside: ", nside) + +# with bench.show("load car kappa"): +# k = enmap.read_map("kappa_car_highres.fits") + +z_min = 0.5 +z_max = 0.51 +mmin = 2e14 +mmax = 2.2e14 + +oras,odecs,ozs,oms = ints.get_agora_halos(z_min = z_min, z_max = z_max, + mass_min = mmin, mass_max = mmax, + massdef='m500') # Get agora halos + +oras = oras[:Nmax] +odecs = odecs[:Nmax] + +print(f"Starting stack on {len(oras)}") +out = 0. +cout = 0. +for i,(ora,odec) in enumerate(zip(oras,odecs)): + coord = np.asarray((odec,ora))*utils.degree + othumb = maps.thumbnail_healpix(kappa,coord,r=rstamp,res=res) # reproject directly from healpix + + # cthumb = reproject.thumbnails(k, coord, r=rstamp, res=res, proj='tan', apod=2*utils.arcmin, + # method="mixed", order=3, oversample=4, pol=None, oshape=None, owcs=None, + # extensive=False, verbose=False, filter=None,pixwin=False,pixwin_order=0) + + + if (i+1)%100==0: print(f"Done {i+1} / {len(oras)}") + out = out + othumb + # cout = cout + cthumb + +N = len(oras) +out = out / N +# cout = cout / N + +lmin = 200 +lmax = 24000 +rmin = 0. +rmax = 10*utils.arcmin +rwidth = 0.5 * utils.arcmin + +cents,kap_agora = lensing.filter_bin_kappa2d(out,lmin=lmin,lmax=lmax,rmin=0.,rmax=rmax,rwidth=rwidth) +# cents,ckap_agora = lensing.filter_bin_kappa2d(cout,lmin=lmin,lmax=lmax,rmin=0.,rmax=rmax,rwidth=rwidth) + + +z = ozs.mean() +ez = ozs.std() +m500c = oms.mean() +em500c = oms.std() + +thetas,kappa_1h,kappa_2h,tot_kappa,cents_t,b1d1h,b1d2h,b1d_t = lensing.kappa_nfw_profiley(mass=m500c,conc=None, + z=z,z_s=1100.,background='critical',delta=500, R_off_Mpc = 0.05, + apply_filter=True, + lmin=lmin,lmax=lmax,res=res,fls=hp.pixwin(nside), # note the healpix pixwin here + rstamp=rstamp,rmin=rmin,rmax=rmax,rwidth=rwidth) +title = f"$M_{{500}}$={m500c/1e14:.2f} $\\pm$ {em500c/1e14:.2f} $10^{{14}} M_{{\\odot}}$ ; z={z:.2f} $\\pm$ {ez:.2f} ; N = {N}" +pl = io.Plotter(xyscale='linlin',xlabel='$\\theta$ (arcmin)',ylabel='$\\kappa$',title=title) +pl.add(cents/utils.arcmin,kap_agora,ls='-',color='r',label=f'Agora stack') +# pl.add(cents/utils.arcmin,ckap_agora,ls='-',color='b',label=f'Agora stack (CAR)') +pl.add(cents_t/utils.arcmin,b1d_t,ls='--',color='k',label='profiley theory') +pl.hline(y=0) +pl.done('profile.png') + +pl = io.Plotter('LCL',xlabel='$\\theta$ (arcmin)',ylabel='$\\theta\\kappa$',title=title) +pl.add(cents/utils.arcmin,kap_agora,ls='-',color='r',label=f'Agora stack') +# pl.add(cents/utils.arcmin,ckap_agora,ls='-',color='b',label=f'Agora stack (CAR)') +pl.add(cents_t/utils.arcmin,b1d_t,ls='--',color='k',label='profiley theory') +pl.hline(y=0) +pl.done('profilet.png') + +pl = io.Plotter(xyscale='linlin',xlabel='$\\theta$ (arcmin)',ylabel='$\\kappa / \\kappa$',title=title) +pl.add(cents/utils.arcmin,kap_agora/b1d_t,ls='-',color='r',label=f'Agora stack / theory') +# pl.add(cents/utils.arcmin,ckap_agora/b1d_t,ls='-',color='b',label=f'Agora stack / theory') +pl.hline(y=1.) +pl._ax.set_ylim(0.3,1.7) +pl.done('profile_rat.png') + + +enmap.write_map("stack.fits",out) \ No newline at end of file From 420c8ef866df58c5cf298fcb7628d214c7eb7830 Mon Sep 17 00:00:00 2001 From: Neha Joshi Date: Wed, 3 Sep 2025 13:17:26 -0400 Subject: [PATCH 10/10] current working files--weird SNR 3 to be checked --- analysis.py | 3 ++- choices.yaml | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/analysis.py b/analysis.py index 82b6473..3d3405f 100644 --- a/analysis.py +++ b/analysis.py @@ -344,7 +344,8 @@ def finish(self, save_transfer=False): pl.add_err(rcents[rsel]/u.arcmin,rmean,yerr=rerr,marker='o') pl.add(self.thetas/u.arcmin, self.kappa) #FIXME: only for flatsky sim? pl._ax.set_xlim(0,10) - pl._ax.set_ylim(-0.01,rmean.max()*1.2) + pl._ax.set_ylim(rmean.min()-0.1,rmean.max()+0.1) + # pl._ax.set_ylim(-0.01,rmean.max()*1.2) pl.hline(y=0) pl.done(self._out('recon_profile.png')) diff --git a/choices.yaml b/choices.yaml index 47dc8bd..38717f4 100644 --- a/choices.yaml +++ b/choices.yaml @@ -46,7 +46,7 @@ ybeam_noise_uk_arcmin: 10.0 skip_lensing: False # This is the profile used for lensing reconstruction in the flatsky sim case -mass: 2.0e+14 +mass: 2.0e+13 delta: 500 #'200' or '500' for critical mass (profiley will calculate m200c or m500c) z: 1.0 theta_max_arc: 20.0 # All templates and reconstructions are masked beyond this radius so that we are not sensitive to the modeling there