In [1]:
import os
os.environ[
"OMP_NUM_THREADS"
] = "64" # for jupyter.nersc.gov otherwise the notebook only uses 2 cores
In [2]:
from pathlib import Path
import healpy as hp
import matplotlib.pyplot as plt
import numpy as np
#import pymaster as nmt
from astropy.io import fits
%matplotlib inline
In [3]:
hp.disable_warnings()
WARNING: AstropyDeprecationWarning: The disable_warnings function is deprecated and may be removed in a future version. [warnings]
In [4]:
plt.style.use("seaborn-talk")
In [5]:
import pysm3 as pysm
import pysm3.units as u
In [6]:
nside = 2048
lmax = 2048
In [7]:
comp = "IQU"
In [8]:
components = list(enumerate(comp))
components
Out[8]:
[(0, 'I'), (1, 'Q'), (2, 'U')]
In [9]:
spectra_components = ["TT", "EE", "BB"]
change this to True if you want to run namaster on notebook
In [10]:
namaster_on_nb = True
In [11]:
datadir=Path("data/")
Setting the inputs¶
Dust maps¶
We use the 2015 GNILC intensity map from the 2nd planck release, as it encodes less contamination from CIB with 21.8' resolution
for Q and U we adopt maps from the 3rd Planck release as they were optimized for polarization studies with 80' reso.
In [12]:
dust_varresI = datadir / "COM_CompMap_Dust-GNILC-F353_2048_21p8acm.fits"
dust_varresP = datadir / "COM_CompMap_IQU-thermaldust-gnilc-varres_2048_R3.00.fits"
In [13]:
if not dust_varresI.exists():
!wget -O $dust_varresI https://portal.nersc.gov/project/cmb/pysm-data/gnilc/inputs/COM_CompMap_Dust-GNILC-F353_2048_21p8acm.fits
In [14]:
if not dust_varresP.exists():
!wget -O $dust_varresP http://pla.esac.esa.int/pla/aio/product-action?MAP.MAP_ID=COM_CompMap_IQU-thermaldust-gnilc-varres_2048_R3.00.fits
Transform maps to double precision for computations
In [15]:
m_planck_varres, h = hp.read_map(dust_varresP, [c + "_STOKES" for c in comp], dtype=np.float64, h=True)
I_planck_varres, h = hp.read_map(dust_varresI, dtype=np.float64, h=True)
Maps from the two releases are in different units MJy/sr the former, and K_CMB the latter, we therefore need to perform some conversion to uK_RJ.
In [16]:
m_planck_varres <<= u.K_CMB
I_planck_varres <<= u.MJy / u.sr
m_planck_varres = m_planck_varres.to(
"uK_RJ", equivalencies=u.cmb_equivalencies(353 * u.GHz)
)
I_planck_varres = I_planck_varres.to(
"uK_RJ", equivalencies=u.cmb_equivalencies(353 * u.GHz)
)
then we are ready to combine both maps into one single TQU map.
In [17]:
m_planck_varres[0] = I_planck_varres
del I_planck_varres
In [18]:
m_planck_varres.dtype
Out[18]:
dtype('float64')
GAL080 Planck mask¶
we mask the galaxy to estimate the power spectra with Namaster before and after small scale injection.
In [19]:
planck_mask_filename = datadir / "HFI_Mask_GalPlane-apo2_2048_R2.00.fits"
if not planck_mask_filename.exists():
!wget -O $planck_mask_filename "http://pla.esac.esa.int/pla/aio/product-action?MAP.MAP_ID=HFI_Mask_GalPlane-apo2_2048_R2.00.fits"
In [20]:
planck_mask = hp.read_map(planck_mask_filename, ["GAL080"])
planck_mask = np.int_(np.ma.masked_not_equal(planck_mask, 0.0).mask)
fsky = planck_mask.sum() / planck_mask.size
print(f"masking {fsky} of the sky")
hp.mollview(planck_mask, title=f"Planck common galactic mask, {comp}")
masking 0.7912631829579672 of the sky
Monopole subtraction¶
Section 2.2 of Planck 2018 XII value reported: 0.13 MJy/sr
we subtract this term only to the I map for the pixels outside the Galactic plane mask.
In [21]:
planck2018_monopole = (0.13 * u.MJy / u.sr).to(
u.uK_RJ, equivalencies=u.cmb_equivalencies(353 * u.GHz)
)
m_planck_varres[0][planck_mask] -= planck2018_monopole
We estimate how many pixels have I< P after we subtract the monopole
In [22]:
maskmono = m_planck_varres[0] ** 2 < m_planck_varres[1] ** 2 + m_planck_varres[2] ** 2
print(
f"{maskmono.sum() } pixels out of { maskmono.size} expected to be NaNs in Log Pol Tens maps "
)
5 pixels out of 50331648 expected to be NaNs in Log Pol Tens maps
In [23]:
plt.figure(figsize=(20, 5))
for i_pol, pol in components:
hp.mollview(
m_planck_varres[i_pol],
title="Planck-GNILC 2058/2018 dust " + pol,
sub=131 + i_pol,
unit=m_planck_varres.unit,
min=-300,
max=300,
)
Transform maps to Poltens quantities¶
In [24]:
def map_to_log_pol_tens(m):
P = np.sqrt(m[1] ** 2 + m[2] ** 2)
log_pol_tens = np.empty_like(m)
log_pol_tens[0] = np.log(m[0] ** 2 - P ** 2) / 2.0
log_pol_tens[1:] = m[1:] / P * np.log((m[0] + P) / (m[0] - P)) / 2.0
return log_pol_tens
def log_pol_tens_to_map(log_pol_tens):
P = np.sqrt(log_pol_tens[1] ** 2 + log_pol_tens[2] ** 2)
m = np.empty_like(log_pol_tens)
exp_i = np.exp(log_pol_tens[0])
m[0] = exp_i * np.cosh(P)
m[1:] = log_pol_tens[1:] / P * exp_i * np.sinh(P)
return m
def sigmoid(x, x0, width, power=4):
"""Sigmoid function given start point and width
Parameters
----------
x : array
input x axis
x0 : float
value of x where the sigmoid starts (not the center)
width : float
width of the transition region in unit of x
power : float
tweak the steepness of the curve
Returns
-------
sigmoid : array
sigmoid, same length of x"""
return 1.0 / (1 + np.exp(-power * (x - x0 - width / 2) / width))
In [25]:
log_pol_tens_varres = map_to_log_pol_tens(m_planck_varres.value)
/tmp/ipykernel_16643/1566809701.py:4: RuntimeWarning: invalid value encountered in log log_pol_tens[0] = np.log(m[0] ** 2 - P ** 2) / 2.0 /tmp/ipykernel_16643/1566809701.py:5: RuntimeWarning: invalid value encountered in log log_pol_tens[1:] = m[1:] / P * np.log((m[0] + P) / (m[0] - P)) / 2.0
check the transformation back and forth
In [26]:
m_back = log_pol_tens_to_map(log_pol_tens_varres)
hp.mollview((m_planck_varres.value - m_back)[1], min=-1e-12, max=1e-12)
del m_back
Checking NaNs on the Poltens map
In [27]:
print(
f"{np.isnan(log_pol_tens_varres[0]).sum() } pixels out of { maskmono.size} are NaNs in Log Pol Tens maps "
)
5 pixels out of 50331648 are NaNs in Log Pol Tens maps
In [28]:
for i in range(3):
log_pol_tens_varres[i, np.isnan(log_pol_tens_varres[i])] = np.nanmedian(
log_pol_tens_varres[i]
)
Set all the NaNs to the map median value
In [29]:
assert np.isnan(log_pol_tens_varres).sum() == 0
In [30]:
hp.write_map(
datadir / "dust_gnilc_logpoltens_varres_nomono.fits",
log_pol_tens_varres,
dtype=np.float32,
overwrite=True,
)
In [31]:
for i_pol, pol in components:
hp.mollview(
log_pol_tens_varres[i_pol],
title="Log Pol tensor " + pol,
sub=131 + i_pol,
unit=m_planck_varres.unit,
)
In [32]:
from scipy.optimize import curve_fit
In [33]:
def model(ell, A, gamma):
out = A * ell ** gamma
out[:2] = 0
return out
In [34]:
def run_anafast(m, lmax):
clanaf = hp.anafast(m, lmax=lmax)
cl = {}
cl["TT"] = clanaf[0]
cl["EE"] = clanaf[1]
cl["BB"] = clanaf[2]
cl["TE"] = clanaf[3]
ell = np.arange(lmax + 1)
cl_norm = ell * (ell + 1) / np.pi / 2
cl_norm[0] = 1
return ell, cl_norm, cl
def run_namaster(m, mask, lmax, nlbins=1):
"""Compute C_ell with NaMaster
Parameters
----------
m : numpy array
T only or TQU HEALPix map
mask : numpy array
mask, 1D, 0 for masked pixels,
needs to have same Nside of the input map
lmax : int
maximum ell of the spherical harmonics transform
Returns
-------
ell : numpy array
array of ell from 0 to lmax (length lmax+1)
cl_norm : numpy array
ell (ell+1)/2pi factor to turn C_ell into D_ell
first element is set to 1
cl : dict of numpy arrays
dictionary of numpy arrays with all components
of the spectra, for now only II, EE, BB, no
cross-spectra
"""
nside = hp.npix2nside(len(mask))
binning = nmt.NmtBin(nside=nside, nlb=nlbins, lmax=lmax, is_Dell=False)
cl = {}
f_0 = nmt.NmtField(mask, [m[0]])
if len(m) == 3:
f_2 = nmt.NmtField(mask, m[1:].copy()) # NaMaster masks the map in-place
cl_namaster = nmt.compute_full_master(f_2, f_2, binning)
cl["EE"] = np.concatenate([[0, 0], cl_namaster[0]])
cl["BB"] = np.concatenate([[0, 0], cl_namaster[3]])
cl_namaster = nmt.compute_full_master(f_0, f_2, binning)
cl["TE"] = np.concatenate([[0, 0], cl_namaster[0]])
elif m.ndim == 1:
m = m.reshape((1, -1))
cl_namaster_I = nmt.compute_full_master(f_0, f_0, binning)
cl["TT"] = np.concatenate([[0, 0], cl_namaster_I[0]])
ell = np.concatenate([[0, 1], binning.get_effective_ells()])
cl_norm = ell * (ell + 1) / np.pi / 2
cl_norm[0] = 1
return ell, cl_norm, cl
In [35]:
print("run anafast on masked sky ")
ell, cl_norm, cl = run_anafast(log_pol_tens_varres, lmax)
run anafast on masked sky
Fitting the power law from the fullsky power spectra¶
We firstly fit a power law spectrum in different multipole ranges. As we get flatter spectral indices for polarization spectra, this will yield to injecting smaller angular scales in TT whose power at given multipole $\ell_*$ is smaller than EE and BB for all $\ell>\ell_*$. To avoid this we therefore force EE and BB small scales to follow the fitted TT power law. We also remark that in this way we get EB ratio closer to ~2.
In [36]:
ell_fit_low = {"TT":100, "EE":30, "BB":30}
ell_fit_high = {"TT":400, "EE":110, "BB":110}
A_fit, gamma_fit, A_fit_std, gamma_fit_std = {},{},{},{}
plt.figure(figsize=(25,5))
for ii, pol in enumerate(spectra_components):
plt.subplot(131+ii)
xdata = np.arange(ell_fit_low[pol], ell_fit_high[pol])
ydata = xdata*(xdata+1)/np.pi/2 * cl[pol][xdata]
(A_fit[pol], gamma_fit[pol]), cov = curve_fit(model, xdata, ydata)
A_fit_std[pol], gamma_fit_std[pol] = np.sqrt(np.diag(cov))
plt.loglog(ell, ell*(ell+1)/np.pi/2 * cl[pol], label="Anafast $C_\ell$")
plt.plot(ell[ell_fit_low[pol]//2:ell_fit_high[pol]*2],
model(ell[ell_fit_low[pol]//2:ell_fit_high[pol]*2], A_fit[pol], gamma_fit[pol]), label="power law fit")
plt.axvline(ell_fit_low[pol], linestyle="--", color="black", label="$ \ell={} $".format(ell_fit_low[pol]))
plt.axvline(ell_fit_high[pol], linestyle="--", color="gray", label="$ \ell={} $".format(ell_fit_high[pol]))
plt.grid()
plt.title(f"{pol} spectrum for dust Dust Pol.Tens " )
plt.ylabel("$\ell(\ell+1)C_\ell/2\pi [\mu K_{RJ}]$")
plt.xlabel(("$\ell$"))
plt.xlim(2, lmax)
print(f"Spectral index from fit for {pol}={gamma_fit[pol]}")
print(f"B-to-E ratio w/ fitted power law at l= {ell_fit_high['BB']} , { A_fit['BB']/A_fit['EE' ]}" )
for ii, pol in enumerate(spectra_components[1:] ):
#we change the EE and BB power laws
A_fit[pol] =A_fit[pol]* ell_fit_high[pol]**( gamma_fit[pol ] - gamma_fit['TT' ] )
gamma_fit[pol] = gamma_fit['TT']
plt.subplot(132+ii)
plt.plot(ell[ell_fit_high[pol] :ell_fit_high[pol]*10],
model(ell[ell_fit_high[pol]:ell_fit_high[pol]*10], A_fit[pol], gamma_fit[pol]),linewidth=3, alpha=.4, color='k',
label="TT power law")
print(f"B-to-E ratio w/ TT power law at l= {ell_fit_high['BB']} , { A_fit['BB']/A_fit['EE' ]}" )
plt.legend()
Spectral index from fit for TT=-1.2872818877318104 Spectral index from fit for EE=-0.33416169211989644 Spectral index from fit for BB=-0.4011846388660788 B-to-E ratio w/ fitted power law at l= 110 , 0.8185188150213931 B-to-E ratio w/ TT power law at l= 110 , 0.5973220042017499
Out[36]:
<matplotlib.legend.Legend at 0x2aad780067c0>
Define Modulation maps¶
as the injected small scales are at different multipoles for intensity and polarization, we consider 2 different modulation maps
Modulation for polarization :¶
- smooth
imap to 5 deg - we saturate all the pixels >4.5 to 1.1
- reduce the dynamic range to range from .1 to 1.1 with MinMax rescaling
Modulation for intensity :¶
- smooth
imap to 5 deg - for the pixels >4.5 MinMax rescaling from 1.1 to 2
- elsewhere MinMax rescaling from .1 to 1.1
In [37]:
ismooth = hp.smoothing(log_pol_tens_varres[0], fwhm=np.radians(5), lmax=lmax)
In [38]:
b1 = 1.1
b2 = 2
a = 0.1
minmax = lambda m, a, b: a + (b - a) * (m - m.min()) / (m.max() - m.min())
modulate_amp_pol = (ismooth) * 1.0
modulate_amp = (ismooth) * 1.0
val = 4.5
mskmd = ismooth > val
modulate_amp_pol[mskmd] = b1
modulate_amp_pol[~mskmd] = minmax(ismooth[~mskmd], a=a, b=b1)
modulate_amp[mskmd] = minmax(ismooth[mskmd], a=b1, b=b2)
modulate_amp[~mskmd] = minmax(ismooth[~mskmd], a=a, b=b1)
In [39]:
del ismooth, mskmd
In [40]:
hp.write_map(datadir / f"modulate_amp_nside{nside}.fits", modulate_amp, dtype=np.float32, overwrite=True)
hp.write_map(
datadir / f"modulate_amp_pol_nside{nside}.fits", modulate_amp_pol, dtype=np.float32, overwrite=True
)
In [41]:
plt.figure(figsize=(15, 5))
hp.mollview(modulate_amp, title="intensity modulation", sub=121)
hp.mollview(modulate_amp_pol, title="polarization modulation", sub=122)
In [42]:
modulate_alm = {}
for name, each_modulate in [("temperature", modulate_amp), ("polarization", modulate_amp_pol)]:
modulate_alm[name] = hp.map2alm(each_modulate, lmax=1.5 * nside, use_pixel_weights=True)
hp.write_alm(datadir / f"gnilc_dust_{name}_modulation_alms_lmax{int(1.5*nside):d}.fits",
modulate_alm[name], overwrite=True, out_dtype=np.float32)
del modulate_amp, modulate_amp_pol, each_modulate
Making maps with small scales¶
- generate the Cl from the fit power spectra and with the sigmoid high-pass filter
- construct Alm from the Cl
- Filter Alm encoding small scales with the sigmoid
- Estimate Alm of the fullsky poltens maps encoding large scales with the sigmoid low-pass filter ,
- Estimate maps from the Alm with large and small scales ,
lsandss - modulate small scales with modulation maps , i.e.
ss' =ss * modulation_maps - define new poltens maps as :
iqu = ls +ss' - transform back to real
IQUmaps
In [43]:
output_nside = 2048 #8192
output_lmax = 2 * output_nside
output_ell = np.arange(output_lmax + 1)
output_cl_norm = output_ell * (output_ell + 1) / np.pi / 2
output_cl_norm[0] = 1
In [44]:
# filter small scales
small_scales_input_cl = [
1
* model(output_ell, A_fit[pol], gamma_fit[pol])
* sigmoid(output_ell, ell_fit_high[pol], ell_fit_high[pol] / 10)
/ output_cl_norm
for pol in spectra_components
]
hp.write_cl(
datadir / f"gnilc_dust_small_scales_logpoltens_cl_lmax{output_lmax}.fits",
small_scales_input_cl,
dtype=np.complex128,
overwrite=True,
)
alm_log_pol_tens_fullsky = hp.map2alm(
log_pol_tens_varres, lmax=lmax, use_pixel_weights=True
)
ii_LS_alm = np.empty_like(alm_log_pol_tens_fullsky)
for ii, pol in enumerate(spectra_components):
sig_func = sigmoid(ell, x0=ell_fit_high[pol], width=ell_fit_high[pol] / 10)
ii_LS_alm[ii] = hp.almxfl(alm_log_pol_tens_fullsky[ii], (1.0 - sig_func) ** 0.2)
np.random.seed(8192)
log_ss = hp.synfast(
small_scales_input_cl + [np.zeros_like(small_scales_input_cl[0])] * 3,
lmax=output_lmax,
new=True,
nside=output_nside,
)
assert np.isnan(log_ss).sum() == 0
modulate_amp = hp.alm2map(modulate_alm["temperature"], output_nside)
log_ss[0] *= modulate_amp
modulate_amp_pol = hp.alm2map(modulate_alm["polarization"], output_nside)
log_ss[1:] *= modulate_amp_pol
assert np.isnan(log_ss).sum() == 0
hp.write_alm(
datadir / f"gnilc_dust_largescale_template_logpoltens_alm_nside{nside}_lmax{int(1.5*nside):d}_complex64.fits",
ii_LS_alm,
lmax=min(int(1.5*nside), output_nside),
out_dtype=np.complex64,
overwrite=True,
)
log_ls = hp.alm2map(ii_LS_alm, nside=output_nside)
ii_map_out = log_ss + log_ls
/tmp/ipykernel_16643/3589521984.py:2: RuntimeWarning: divide by zero encountered in power out = A * ell ** gamma
In [45]:
del log_ls, ii_LS_alm
In [46]:
output_map = log_pol_tens_to_map(ii_map_out)
hp.write_map(
datadir / f"gnilc_dust_template_nside{output_nside}_float32.fits",
output_map,
dtype=np.float32,
overwrite=True,
)
In [47]:
output_map_alm = hp.map2alm(output_map, lmax=output_lmax)
In [ ]:
hp.write_alm(
datadir / f"dust_gnilc_template_alm_nside{output_nside}_lmax{output_lmax}_complex64.fits",
output_map_alm,
out_dtype = np.complex64,
overwrite=True,
)
Plot maps¶
In [49]:
lat=15
plt.figure(figsize=(15,10))
hp.gnomview(ii_map_out[0] ,cmap='RdBu', title='i w/ small scales ', rot=[0,lat],reso=3.75,xsize=320, sub=234 )
hp.gnomview(log_pol_tens_varres [0],cmap='RdBu', title='i', rot=[0,lat],reso=3.75,xsize=320, sub=231 )
hp.gnomview((modulate_amp),cmap='RdBu', title=' modulation I ', rot=[0,lat],reso=3.75,xsize=320, sub=233, )
hp.gnomview((m_planck_varres [0]),cmap='RdBu', title=' I GNILC ', rot=[0,lat],reso=3.75,xsize=320, sub=232, )
hp.gnomview((log_ss)[0],cmap='RdBu', title=' small scales ', rot=[0,lat], reso=3.75,xsize=320, sub=236, )
hp.gnomview(output_map[0],cmap='RdBu' , title='I w/ small scales ', rot=[0,lat],reso=3.75,xsize=320, sub=235 )
plt.figure(figsize=(15,10))
hp.gnomview(ii_map_out[1],cmap='RdBu' , title='q w/ small scales ', rot=[0,lat],reso=3.75,xsize=320, sub=234 )
hp.gnomview(log_pol_tens_varres [1],cmap='RdBu', title='q', rot=[0,lat],reso=3.75,xsize=320, sub=231 )
hp.gnomview((modulate_amp_pol),cmap='RdBu', title=' modulation ', rot=[0,lat],reso=3.75,xsize=320, sub=233, )
hp.gnomview((m_planck_varres [1]),cmap='RdBu', title=' Q GNILC ', rot=[0,lat],reso=3.75,xsize=320, sub=232, )
hp.gnomview((log_ss)[1],cmap='RdBu', title=' small scales ', rot=[0,lat], reso=3.75,xsize=320, sub=236, )
hp.gnomview(output_map[1] ,cmap='RdBu', title='Q w/ small scales ', rot=[0,lat],reso=3.75,xsize=320, sub=235 )
In [50]:
del log_ss, ii_map_out, modulate_amp, modulate_amp_pol
In [51]:
spectra_components += ["TE"]
In [52]:
output_planck_mask = hp.ud_grade(planck_mask, output_nside)
In [53]:
hp.mollview(output_planck_mask)
Validate outputs estimating NaMaster power spectra on both input and outputs¶
we have expanded our validation by means of 3 figures of merit:
- TT,EE, BB power spectra with Namaster
- B-to-E ratio
- Polarization fraction
Moreover, we have considered 4 masks to evaluate the quality of the small scales injected for different $f_{sky} = 0.02,0.2, 0.4,0.8 $.
Another small difference with previous analysis is that we estimated the spectra on binned equally-spaced multipoles. We choose $\Delta \ell =35, 25,15,4$ respectively for the 4 masks.
In [54]:
bk15_mask_filename = datadir/ "BK15_region_Gal_apo.fits"
if not bk15_mask_filename.exists():
!wget -O $bk15_mask_filename https://portal.nersc.gov/project/cmb/pysm-data/gnilc/inputs/BK15_region_Gal_apo.fits
In [55]:
planck_masks = {k.lower():hp.read_map(datadir/ "HFI_Mask_GalPlane-apo2_2048_R2.00.fits", [k]) for k in ["GAL{:03d}".format(frac) for frac in [20,40,80]]}
planck_masks["BK"] = hp.read_map(bk15_mask_filename)
planck_masks ={ k:np.ma.masked_equal(m , 1 ).mask for k,m in planck_masks.items() }
In [56]:
nlb= {'BK':35, 'gal020':25, 'gal040':15,'gal080':4 }
In [57]:
import pymaster as nmt
In [58]:
plt.figure(figsize=(10,5))
for jj, k in enumerate(planck_masks.keys()) :
fspectra = datadir / f"dust_gnilc_hybrid_out_nside2048_float32_{k}_spectra.npz"
if os.path.exists(fspectra):
print("read Namaster spectra ")
output_ell = np.load(fspectra)["ell"]
cl_out = {kk: np.load(fspectra)[kk] for kk in spectra_components}
elif namaster_on_nb:
print("run Namaster ")
output_ell, output_cl_norm, cl_out = run_namaster(
output_map, mask=output_planck_mask, lmax=output_lmax, nlbins = nlb[k]
)
np.savez(fspectra, ell=output_ell, cl_norm=output_cl_norm, **cl_out)
#else:
# print("run anafast on masked sky ")
# output_ell, cl_norm, cl_out = run_anafast(output_map * planck_mask, lmax)
fspectra = datadir / f"dust_gnilc_varres_no_monopole_{k}_spectra.npz"
if os.path.exists(fspectra):
print("read Namaster spectra ")
input_ell = np.load(fspectra)["ell"]
cl_in = {kk: np.load(fspectra)[kk] for kk in spectra_components}
elif namaster_on_nb:
print("run Namaster ")
input_ell, input_cl_norm, cl_in = run_namaster(
m_planck_varres, mask=planck_mask, lmax=output_lmax, nlbins = nlb[k]
)
np.savez(fspectra, ell=input_ell, cl_norm=input_cl_norm, **cl_in)
#else:
# print("run anafast on masked sky ")
# ell, cl_norm, cl_out = run_anafast(output_map * planck_mask, lmax)
plt.subplot(2,2,jj+1 )
for ii, pol in enumerate(["TT","EE", "BB" ]):
if jj==3:
plt.loglog(output_ell, cl_out [pol], color='C%d'%ii, label=pol
)
else:
plt.loglog(output_ell, cl_out [pol], color='C%d'%ii, )
plt.loglog(input_ell, cl_in [pol], color='C%d'%ii,alpha=.5,)
plt.grid()
plt.legend(title=(k +" mask "), fontsize=15)
plt.ylabel("$ C_\ell [\mu K_{RJ}]$")
plt.xlabel(("$\ell$"))
plt.xlim(10,2e3)
plt.tight_layout()
run Namaster run Namaster
No handles with labels found to put in legend. No handles with labels found to put in legend. No handles with labels found to put in legend.
run Namaster
No handles with labels found to put in legend. No handles with labels found to put in legend. No handles with labels found to put in legend.
run Namaster run Namaster
No handles with labels found to put in legend. No handles with labels found to put in legend. No handles with labels found to put in legend.
run Namaster run Namaster
In [59]:
plt.figure(figsize=(10,5))
for jj, k in enumerate(planck_masks.keys()) :
fspectra = datadir / f"dust_gnilc_hybrid_out_nside2048_float32_{k}_spectra.npz"
if os.path.exists(fspectra):
print("read Namaster spectra ")
output_ell = np.load(fspectra)["ell"]
cl_out = {kk: np.load(fspectra)[kk] for kk in spectra_components}
elif namaster_on_nb:
print("run Namaster ")
output_ell, output_cl_norm, cl_out = run_namaster(
output_map, mask=output_planck_mask, lmax=output_lmax, nlbins = nlb[k]
)
np.savez(fspectra, ell=output_ell, cl_norm=output_cl_norm, **cl_out)
#else:
# print("run anafast on masked sky ")
# output_ell, cl_norm, cl_out = run_anafast(output_map * planck_mask, lmax)
fspectra = datadir / f"dust_gnilc_varres_no_monopole_{k}_spectra.npz"
if os.path.exists(fspectra):
print("read Namaster spectra ")
input_ell = np.load(fspectra)["ell"]
cl_in = {kk: np.load(fspectra)[kk] for kk in spectra_components}
elif namaster_on_nb:
print("run Namaster ")
input_ell, input_cl_norm, cl_in = run_namaster(
m_planck_varres, mask=planck_mask, lmax=output_lmax, nlbins = nlb[k]
)
np.savez(fspectra, ell=input_ell, cl_norm=input_cl_norm, **cl_in)
#else:
# print("run anafast on masked sky ")
# ell, cl_norm, cl_out = run_anafast(output_map * planck_mask, lmax)
plt.subplot(2,2,jj+1 )
plt.semilogx(output_ell, cl_out ["BB"]/cl_out["EE"] , color='C%d'%0, )
plt.semilogx(input_ell, cl_in ["BB"]/cl_in["EE"], color='C%d'%0,alpha=.5,)
plt.grid()
plt.legend(title=(k +" mask "), fontsize=15)
plt.ylabel(("B-to-E ratio"))
plt.xlabel(("$\ell$"))
plt.ylim(0,1)
plt.xlim(10,2e3 )
plt.tight_layout()
/tmp/ipykernel_16643/2708351735.py:35: RuntimeWarning: invalid value encountered in true_divide plt.semilogx(output_ell, cl_out ["BB"]/cl_out["EE"] , color='C%d'%0, ) /tmp/ipykernel_16643/2708351735.py:36: RuntimeWarning: invalid value encountered in true_divide plt.semilogx(input_ell, cl_in ["BB"]/cl_in["EE"], color='C%d'%0,alpha=.5,) No handles with labels found to put in legend. /tmp/ipykernel_16643/2708351735.py:35: RuntimeWarning: invalid value encountered in true_divide plt.semilogx(output_ell, cl_out ["BB"]/cl_out["EE"] , color='C%d'%0, ) /tmp/ipykernel_16643/2708351735.py:36: RuntimeWarning: invalid value encountered in true_divide plt.semilogx(input_ell, cl_in ["BB"]/cl_in["EE"], color='C%d'%0,alpha=.5,) No handles with labels found to put in legend. /tmp/ipykernel_16643/2708351735.py:35: RuntimeWarning: invalid value encountered in true_divide plt.semilogx(output_ell, cl_out ["BB"]/cl_out["EE"] , color='C%d'%0, ) /tmp/ipykernel_16643/2708351735.py:36: RuntimeWarning: invalid value encountered in true_divide plt.semilogx(input_ell, cl_in ["BB"]/cl_in["EE"], color='C%d'%0,alpha=.5,) No handles with labels found to put in legend. /tmp/ipykernel_16643/2708351735.py:35: RuntimeWarning: invalid value encountered in true_divide plt.semilogx(output_ell, cl_out ["BB"]/cl_out["EE"] , color='C%d'%0, ) /tmp/ipykernel_16643/2708351735.py:36: RuntimeWarning: invalid value encountered in true_divide plt.semilogx(input_ell, cl_in ["BB"]/cl_in["EE"], color='C%d'%0,alpha=.5,) No handles with labels found to put in legend.
read Namaster spectra read Namaster spectra read Namaster spectra read Namaster spectra read Namaster spectra read Namaster spectra read Namaster spectra read Namaster spectra
Estimate Pol. fraction distribution¶
In [60]:
get_polfrac = lambda m: np.sqrt(m[1] ** 2 + m[2] ** 2) / m[0]
In [61]:
Pout = get_polfrac(output_map)
Pin = get_polfrac(m_planck_varres.value)
In [62]:
plt.figure(figsize=(15, 5))
hp.mollview(
np.log10(Pin), title=" Pol.Frac input", sub=121, min=-3, max=0, unit="Log10(p )"
)
hp.mollview(
np.log10(Pout), title="Pol.Frac output", sub=122, min=-3, max=0, unit="Log10(p)"
)
In [63]:
logpin = np.log10(Pin)
logpout = np.log10(Pout)
In [64]:
plt.figure(figsize=(10,5))
for jj,pm in enumerate(planck_masks.items()) :
k =pm[0]
msk =pm[1]
plt.subplot(2,2,jj+1 )
h,edg= np.histogram(
logpout[msk] ,
bins=np.linspace(-4, 0, 100), density=True)
xb = np.array([(edg[i] +edg[i+1])/2 for i in range(edg.size-1)])
plt.plot( xb,h,
lw=3,
color='C0', alpha=.5,label='output'
)
h, edg= np.histogram(
logpin[msk] ,
density=True ,
bins=np.linspace(-4, 0, 100))
xb = np.array([(edg[i] +edg[i+1])/2 for i in range(edg.size-1)])
plt.plot(xb,h ,
lw=3,alpha=.5 ,color='C0' ,linestyle='--',label='input'
)
plt.ylabel("norm.counts", fontsize=14)
plt.xlabel(r"$\log10( p )$", fontsize=14)
plt.legend(title=(k +" mask "), fontsize=13, loc='best' )
plt.xlim(-3,-0.5)
plt.tight_layout()
