zonca commented 2 years ago

The implementation of the GNILC-based dust templates is complete, see details and the executed notebook at:

https://github.com/galsci/pysm/pull/97#issuecomment-1006313372

Now it would be useful to have some independent validation, load the maps, take spectra, compare known regions, check everything is reasonable (see also #100)

zonca commented 2 years ago

Files location

https://portal.nersc.gov/project/cmb/pysm-dev/gnilc/

or at NERSC in:

/global/project/projectdirs/cmb/www/pysm-dev/gnilc

zonca commented 2 years ago

Available files

Output templates

In uK_RJ (see also in FITS header)

gnilc_dust_template_nside2048_float32.fits
gnilc_dust_template_nside4096_float32.fits
gnilc_dust_template_nside8192_float32.fits

Input GNILC products

For temperature:

COM_CompMap_Dust-GNILC-F353_2048_21p8acm.fits

For polarization:

COM_CompMap_IQU-thermaldust-gnilc-varres_2048_R3.00.fits

Masks

BK15_region_Gal_apo.fits
HFI_Mask_GalPlane-apo2_2048_R2.00.fits

Intermediate products

Large scale Alms in logpoltens

gnilc_dust_largescale_template_logpoltens_alm_nside2048_lmax3072_complex64.fits

Small scales power spectrum

gnilc_dust_small_scales_logpoltens_cl_lmax16384.fits

Modulation of the small scales

gnilc_dust_temperature_modulation_alms_lmax3072.fits
gnilc_dust_polarization_modulation_alms_lmax3072.fits

zonca commented 2 years ago

@brandonshensley @seclark @NicolettaK @yaojian95 @giuspugl @delabru

zonca commented 2 years ago

Spectral index and Dust temperature

At the same location above I have also copied the outputs of the GNILC beta and Td processing, see the notebooks in https://github.com/galsci/pysm/pull/104 for more details.

Output maps

6 maps, 2 componets at 3 resolutions:

gnilc_dust_[beta,Td]_nside[2048,4096,8192]_float32.fits

Intermediate products

gnilc_dust_largescale_template_beta_alm_nside2048_lmax1024_complex64.fits
gnilc_dust_largescale_template_Td_alm_nside2048_lmax1024_complex64.fits
gnilc_dust_small_scales_beta_cl_lmax16384.fits
gnilc_dust_small_scales_Td_cl_lmax16384.fits

For modulation I reuse the "temperature" modulation above.

b-thorne commented 2 years ago

A while ago I was calculating power spectra of the new log pol tens templates. I found that the power spectra of smaller patches of sky exhibits more of a break at the transition scale around ell of 100 on smaller patches of sky. For example, the bottom panel of this plot:

Screen Shot 2022-03-17 at 11 52 59 AM

The bottom panel shows BB, and the green line shows the spectra calculated on the BICEP / Keck patch, with the dashed line being the input map COM_CompMap_Dust-GNILC-F353_2048_21p8acm.fits, and the solid line the derived log pol tens template, at the time. Initially I was concerned about this plot because we would be under predicting dust BB in the BK patch in the important ell range of 100-300 by almost an order of magnitude.

Andrea and Giuseppe were unable to reproduce these plots, and so I reran the spectra, using the updated log pol tens templates in this thread, as well as the BK mask in this thread: BK15_region_Gal_apo.fits. This is a comparison of the various masks (Andrea's BK mask is labelled BK_NERSC):

validation_masks (1)

With the updated templates, I find that the discrepancy in the BK patch (both my original mask, and the new mask), is reduced by a factor of a few, whilst the BB spectrum in the LR15 region stays a little low (and strangely lower than the two BK regions). I've also plotted in the BB panel the maximum likelihood BK dust model for comparison.

spectra_by_field_ana_NSIDE1024 (7)

Based on this final plot, I'm not sure this is an issue worth investigating much further, but I wanted to record it here.

giuspugl commented 2 years ago

Thanks a lot @bthorne93 for this follow up! it is great to see this validation with the latest updated models!

Firstly a comment on what has changed between the two versions of plots is that in the former one we had injected small scales using the fitted power law indices from TT, EE and BB spectra (respectively about -1.2, -0.3,-0.4 ). In the latter, the small scales are included adopting the same TT index for all the IQU maps. This results in bumping up EE and BB modes in the final spectra (so the edge effect is milder than the former plot) , but it is not intuitive to me how this translate into effectively having less power in LR15 than BK patches . Any clue?
about the edge we observe especially in EE and BB from LR34 and below i think that we can't do much more than what we currently have done.. Although it 's great to see that from the spectra going down very well with a power law (remember we modulate the small scales so that we have something reliable at fsky=80%) up to fsky=20%, i was expecting some edge to become more visible for smaller patches . This points out that we really need to implement a modulation methodology that encodes also local info ( with needlets or something similar) .
lastly thanks a lot for including the solid black line from BK maximum likelihood. It adds an extra-data point in our understanding: the BK estimates of dust were essentially driven by the dust spectrum from GNILC maps at ell<100, the spectral index seems to me kind of parallel to the dashed green and orange lines. However, the BK estimated amplitude seems to me factor of 2 lower than the one derived from the GNILC template. Do you know why ?

b-thorne commented 2 years ago

@giuspugl

I am not sure how fixing the power law index to the TT value, which is steeper than the other two values you showed, would bump up the EE/BB spectra at multipoles higher than the pivot scale (which I think is ~100, right?). I would have thought that change would actually make the simulated power lower at scales above the pivot scale.
Since we are fixing the slope of the power spectrum to the measured TT value of -1.2, I think it makes sense that we see a much steeper slope in BB than BICEP / Keck measured. The BICEP / Keck dataset has higher S/N on the dust amplitude in their 220 GHz channel, rather than in the Planck 353 GHz data. Perhaps in this region we have low S/N, resulting in slightly biased autospectra. This may also explain why there is very little difference between the BK and LR15 regions?

brandonshensley commented 2 years ago

@giuspugl The issue of alpha_BB is being raised, so I think it may be worth digging back into @bthorne93's question 1 above.

giuspugl commented 2 years ago

I am not sure how fixing the power law index to the TT value, which is steeper than the other two values you showed, would bump up the EE/BB spectra at multipoles higher than the pivot scale (which I think is ~100, right?). I would have thought that change would actually make the simulated power lower at scales above the pivot scale.

@bthorne93 , sorry I totally lost track of this comment! I totally agree that choosing a steeper spectr. index would lower the power than a flat one. (i think there was a typo in my previous comment). A solution we can discuss is to set 2 pivot scales (ell_p1 ~100 and ell_p2 ~600 ) :

if ell<ell_p1 : we use GNILC templates,
if ell_p1<ell <ell_p2 use flat spectral index consistent w/ BK and Planck data
if ell> ell_p2 use a steep spectral index to avoid crossing of EE and BB over TT

This approach is still in line with our philosophy, i.e. use as much as we could observations and extend simulations with physical assumptions.

seclark commented 2 years ago

I like this approach. The steeper spectral index should be chosen to anticipate high nside. Maybe we select the index to ensure that EE, BB < TT out to 8192 at least.

giuspugl commented 2 years ago

In this notebook you'll find an attempt to implement what i proposed above:

if ell<ell_p1 : we use GNILC templates,

if ell_p1<ell <ell_p2 use flat spectral index consistent w/ BK and Planck data

if ell> ell_p2 use a steep spectral index to avoid crossing of EE and BB over TT

Even in this case, our benchmark is to match small and large scales for the mask GAL080, however this

/!\ This new implementation employs several new features wrt the previous d9, d10 models:

two pivot scales ell_p1,ell_p2 for small scale injection
spectral indices for EE, BB coming from literature , Planck 2018 XI and for TT from Miville-Deschenes 2016
we inject small scales with non-zero TE from Planck 2018 XI
Inside the GAL097 mask (i.e. along the Gal. midplane) we don't inject small scales, we simply keep the ones observed at high SNR by Planck.

Spectra at intermediate and low Gal. latitudes :

Spectra at high Gal. latitudes

giuspugl commented 2 years ago

Updates on the models

I 've finalized a new model of dust with small scales employing several changes wrt what has been done before:

two pivot scales ell1=110,ell2=800 for small scale injection (see above)
spectral indices for EE, BB coming from literature , Planck 2018 XI and for TT from Miville-Deschenes 2016
we inject small scales with non-zero TE from Planck 2018 XI
Inside the GAL097 mask (i.e. along the Gal. midplane) we don't inject small scales, we simply keep the ones observed at high SNR by Planck.
no modulation is employed for polarization small scales. the modulation is applied only in the T map. As the injection is performed in the log pol. tensor domain, fluctuations in the total polarization amplitude are modulated by the ones in T map. the remaining fluctuations due to turbulence in the magnetic field, are roughly isotropic on the sky so don’t need to be modulated.

Power Spectrum indices:

Maps

New maps looks noisier compare to the templates or the ones previously produced in pysm3, however this is kind of expected, small scales contribute more in Q and U maps when a flatter spectral index is employed.

Power Spectra in `HFI Masks`

More importantly, avoiding modulation in P actually helps in improving the quality of the spectra estimated across several masks:

Power Spectra in South-pole patches:

Polarization fraction

Finally, we do observe a mild increase around the peak of polarization fraction distribution, but the small scale injection do not alter the overall shape of the distribution .

Notebook and maps

Notebook is publicly available : https://gist.github.com/giuspugl/3f63453d9d6351bdd4da64f6560fb934 New maps are available in a shared folder at NERSC: /global/cscratch1/sd/giuspugl/cmbs4/pysm_validation/

brandonshensley commented 2 years ago

Thanks @giuspugl, this is great! One question though: why does the new T map look like it retains less information from the input map relative to the old T map?

giuspugl commented 2 years ago

why does the new T map look like it retains less information from the input map relative to the old T map?

My understanding is that small scale content in P map is larger especially at high multipoles (because it 's less modulated and because E and B mode spectral indices are flatter) therefore you get a larger loss of info also in the previous T map. Remember that those small scales are mixed when you transform back from poltens to real IQU quantities . In particular I = e^i cosh p.

zonca commented 2 years ago

@brandonshensley can someone do a further validation with @giuspugl's maps at /global/cscratch1/sd/giuspugl/cmbs4/pysm_validation/ before I implement this into PySM?

brandonshensley commented 2 years ago

Yes, let's please hold off for now in implementation in PySM. Thanks!

seclark commented 2 years ago

Agree on holding off for now. Are there validation measurements needed that don't yet have someone assigned?

zonca commented 2 years ago

@giuspugl apologies, I lost the last 10 min of today's call. Are you going to post the notebook here?

giuspugl commented 2 years ago

Updates on dust

Below a summary of the updates implemented on new dust models :

input templates from GNILC variable resolution and unires GNILC maps
two pivot scales ell1=110,ell2=2000 for small scale injection (see my previous posts )
spectral indices for EE, BB coming from literature , Planck 2018 XI and for TT from Miville-Deschenes 2016
we inject small scales with non-zero TE from Planck 2018 XI
Inside the GAL097 mask (i.e. along the Gal. midplane) we don't inject small scales, we simply keep the ones observed at high SNR by Planck.
Modulation of qu maps with a single map, p: i. to avoid modulation w/ negative values, ii. to preserve non-zero TE
we propose to modulate small scales as it has been done in pysm2(https://arxiv.org/pdf/1608.02841.pdf) , with a couple of differences: i. small scales are expected to be injected with non-gaussian content (thanks to the logpoltens formalism); ii. split the sky with high reso pixels (nside=8) ; iii. use amplitude of E-mode spectra to derive the modulation template

Validation

TT, TE, EE, BB Power Spectra at different GAL masks

E-to- B ratio at different GAL masks

EE BB spectra in small patches

Polarization Fraction pixel distribution

Notebook and final dust maps

you can find my final notebook here Maps : /global/cscratch1/sd/giuspugl/cmbs4/pysm_validation/dust_gnilc_varres_pysmmod.fits

zonca commented 2 years ago

thanks @giuspugl! I'll start implementing this into PySM

zonca commented 2 years ago

@giuspugl can you please make the data folder /global/cscratch1/sd/giuspugl/workstation/FG_extension/extending_gnilc_dust/ readable?

zonca commented 2 years ago

@giuspugl can you please make the data folder /global/cscratch1/sd/giuspugl/workstation/FG_extension/extending_gnilc_dust/ readable?

nevermind, I reran namaster on the patches

zonca commented 2 years ago

@giuspugl clarification on cell 45:

    ii_LS_alm[ii] = hp.almxfl(alm_log_pol_tens_fullsky[ii], (1.0 - sig_func) ** 0.2)

if we are modulating the small scales in C_ell with sig_func, shouldn't we modulate these with sqrt(1 - sig_func)?

giuspugl commented 2 years ago

@giuspugl clarification on cell 45:
    ii_LS_alm[ii] = hp.almxfl(alm_log_pol_tens_fullsky[ii], (1.0 - sig_func) ** 0.2)
if we are modulating the small scales in C_ell with sig_func, shouldn't we modulate these with sqrt(1 - sig_func)?

This is something we found in the past to have a rounder transition at around ell1 in the sigmoid, unfortunately the sqrt was too sharp.

zonca commented 2 years ago

ok, I have modified giuspugl notebook to have a fixed seed and used it to compare with my 2-step implementation (first generate inputs in spherical harmonics space, basically inputs to d11, then from that create the d10 templates), for reference:

https://gist.github.com/c5f71634f4466a52945eb22b0044a3c2 https://gist.github.com/3bf14002bb54d2959c9cb02de48674fe

zonca commented 2 years ago

ok, implemented in #133

galsci / pysm

Validation of GNILC dust templates #101

Files location

Available files

Output templates

Input GNILC products

Masks

Intermediate products

Spectral index and Dust temperature

Output maps

Intermediate products

Spectra at intermediate and low Gal. latitudes :

Spectra at high Gal. latitudes

Updates on the models

Power Spectrum indices:

Maps

Power Spectra in `HFI Masks`

Power Spectra in South-pole patches:

Polarization fraction

Notebook and maps

Updates on dust

Validation

TT, TE, EE, BB Power Spectra at different GAL masks

E-to- B ratio at different GAL masks

EE BB spectra in small patches

Polarization Fraction pixel distribution

Notebook and final dust maps

galsci / pysm

Validation of GNILC dust templates #101

Files location

Available files

Output templates

Input GNILC products

Masks

Intermediate products

Spectral index and Dust temperature

Output maps

Intermediate products

Spectra at intermediate and low Gal. latitudes :

Spectra at high Gal. latitudes

Updates on the models

Power Spectrum indices:

Maps

Power Spectra in HFI Masks

Power Spectra in South-pole patches:

Polarization fraction

Notebook and maps

Updates on dust

Validation

TT, TE, EE, BB Power Spectra at different GAL masks

E-to- B ratio at different GAL masks

EE BB spectra in small patches

Polarization Fraction pixel distribution

Notebook and final dust maps

Power Spectra in `HFI Masks`