Debugging gaussian scales for GNILC dust

View / edit / reply to this conversation on ReviewNB

zonca commented on 2021-05-11T00:59:48Z ----------------------------------------------------------------

Jacques suggested to use the apodized mask, I was using the 2 degrees mask from planck in the last version of my notebook https://nbviewer.jupyter.org/gist/zonca/b398f85f701088ba60e3e66a20e0f48f.

Unless we are reason to believe the non apodized mask is better.

giuspugl commented on 2021-05-11T21:31:08Z ----------------------------------------------------------------

Ok will do that, though not expecting large differences from this .

review-notebook-app[bot] commented 3 years ago

View / edit / reply to this conversation on ReviewNB

zonca commented on 2021-05-11T00:59:49Z ----------------------------------------------------------------

as we discussed, another option is to use the unires version for fitting, the worry is that the variable resolution could have extra signal at small scales coming from some regions of the sky. I think it would be nice to have 2 versions of the notebook to compare the spectra.

giuspugl commented on 2021-05-11T21:32:03Z ----------------------------------------------------------------

Yeah, I did also that case, will include a new notebook with the unires version

review-notebook-app[bot] commented 3 years ago

View / edit / reply to this conversation on ReviewNB

zonca commented on 2021-05-11T00:59:49Z ----------------------------------------------------------------

what about just removing the value from the Planck paper? does it give NaNs?

It probably makes more sense to adjust the monopole before transforming to poltens.

giuspugl commented on 2021-05-11T22:46:37Z ----------------------------------------------------------------

If we subtract the value reported by Planck the NaNs reduce to 149 pixels, think we should be fine

review-notebook-app[bot] commented 3 years ago

View / edit / reply to this conversation on ReviewNB

zonca commented on 2021-05-11T00:59:50Z ----------------------------------------------------------------

I think it is useful to compare with anafast but we should deconvolve the mask with polspice or namaster

giuspugl commented on 2021-05-11T21:36:32Z ----------------------------------------------------------------

I agree with you, the results shown in the plots are obtained with namaster, sorry for the confusion

review-notebook-app[bot] commented 3 years ago

View / edit / reply to this conversation on ReviewNB

zonca commented on 2021-05-11T00:59:51Z ----------------------------------------------------------------

how are those spectra computed? by a previous run?

it would be useful to make an if statement. If the npz is available, load it, otherwise compute it and save it.

giuspugl commented on 2021-05-11T21:37:41Z ----------------------------------------------------------------

they are computed with namaster, but unfortunately i can't run namaster on jupyter notebook , so have to run it on an interactive node .

review-notebook-app[bot] commented 3 years ago

View / edit / reply to this conversation on ReviewNB

zonca commented on 2021-05-11T00:59:51Z ----------------------------------------------------------------

In the last version of #72 I reduced those ell ranges, I think they are too high here, we are doing the fit where the effect of the beam is very large.

giuspugl commented on 2021-05-11T21:40:07Z ----------------------------------------------------------------

ok will perform a fit at lower ells

zonca commented 3 years ago

thanks Giuseppe! I added a few comments, but better also double-check my suggestions with someone more expert in this kind of analysis!

giuspugl commented 3 years ago

Ok will do that, though not expecting large differences from this .

View entire conversation on ReviewNB

giuspugl commented 3 years ago

Yeah, I did also that case, will include a new notebook with the unires version

View entire conversation on ReviewNB

giuspugl commented 3 years ago

I agree with you, the results shown in the plots are obtained with namaster, sorry for the confusion

View entire conversation on ReviewNB

giuspugl commented 3 years ago

they are computed with namaster, but unfortunately i can't run namaster on jupyter notebook , so have to run it on an interactive node .

View entire conversation on ReviewNB

giuspugl commented 3 years ago

ok will perform a fit at lower ells

View entire conversation on ReviewNB

giuspugl commented 3 years ago

If we subtract the value reported by Planck we don't have the NaNs . So will use that value and remove the monopole before hand

View entire conversation on ReviewNB

brandonshensley commented 3 years ago

The argument against the unires maps is that treatment of polarization systematics in Planck evolved from the 2015 to the 2018 release. I wonder if that accounts for the qualitative differences we are seeing? Rather than use the unires QU maps, it may be worth investigating the NPIPE raw 353 GHz QU maps instead.

giuspugl commented 3 years ago

I addressed most of @zonca comments on the notebooks, in particular :

subtracting the monopole term in the intensity map (not in the log pol tens transformed quantities) with the value reported in Planck 2018
performing the whole analysis with the varres maps with the GAL080 mask apodized with 2 degree radius
namaster spectra
fit the power law at the ranges as in #72 , i.e. for TT ell =400-600, for EE and BB ell=30-110 I am not finding flat spectral indices , see the spectra below
The small scales are coadded in the logpol tens maps and the namaster spectra are shown below

notice in particular the light green line ( with large scale only) and the orange line encoding spectra with gaussian small scales injected . Finally, I transformed back the maps into real quantities and got reasonable power spectra (mild decrease in power in TT at ell>600)

gonna clean up a bit the notebook and commit it to this PR.

giuspugl commented 3 years ago

The argument against the unires maps is that treatment of polarization systematics in Planck evolved from the 2015 to the 2018 release. I wonder if that accounts for the qualitative differences we are seeing? Rather than use the unires QU maps, it may be worth investigating the NPIPE raw 353 GHz QU maps instead.

discussed with @keskitalo and we agreed that this is something that can be surely done, although we have to deal with synchrotron and CMB residuals when using the NPIPE maps as those haven't been processed with component separation.

brandonshensley commented 3 years ago

If everything is working for varres--and I think it looks really good!--no need to go to the NPIPE maps. Relieved to see those -0.5 spectra!

giuspugl commented 3 years ago

Notebook with varres maps available at https://gist.github.com/giuspugl/4fb30306e99dfed4a193c80fe9c3c94a

review-notebook-app[bot] commented 3 years ago

View / edit / reply to this conversation on ReviewNB

zonca commented on 2021-05-13T18:39:05Z ----------------------------------------------------------------

I think we want to keep the mask apodized, that should provide better behavior in namaster

giuspugl commented on 2021-05-13T18:48:09Z ----------------------------------------------------------------

Yes sorry if that wasn't clear, when running namaster always the mask is apodized.

review-notebook-app[bot] commented 3 years ago

View / edit / reply to this conversation on ReviewNB

zonca commented on 2021-05-13T18:39:06Z ----------------------------------------------------------------

Just for temperature, we should probably multiply the small scales by a map of the intensity heavily smoothed, maybe 3 degrees? And normalized, I think we need to adjust the normalization so that the spectra look reasonable (have no sharp boundaries or other anomalous features) both when we take spectra on the masked sky and on the full sky.

review-notebook-app[bot] commented 3 years ago

View / edit / reply to this conversation on ReviewNB

zonca commented on 2021-05-13T18:39:07Z ----------------------------------------------------------------

it is really strange this big loss of power, is it possible that the small scales were lost in the processing?

Maybe overplotting just the small scales both in logpoltens and TT can help understand what is going on.

zonca commented 3 years ago

thanks @giuspugl! great progress, I am worried that the loss of power in TT (IQU) compared to Poltens could be caused by a bug somewhere.

giuspugl commented 3 years ago

Yes sorry if that wasn't clear, when running namaster always the mask is apodized.

View entire conversation on ReviewNB

seclark commented 3 years ago

Agreed, looking great @giuspugl ! Re: loss of power in TT: is there also a small knee in BB at the same ell?

giuspugl commented 3 years ago

I've dug a bit into the mismatch we discussed last week, see the gist https://gist.github.com/giuspugl/4fb30306e99dfed4a193c80fe9c3c94a I am finding that the reason for the mismatch is mostly driven by the fact that the i map with large scales in the logpol tensor formalism shows the largest contribution when we transform back to IQU maps. I think we are hitting a limitation of this proposed formalism proposed, although i 'd like to stress that a by-product of this formalism is that the smaller scales even if simulated uniformly they present a modulation by the larger scales of the intensity maps.

brandonshensley commented 3 years ago

If we simply let i = ln(I) does everything look ok, or are there other problems that fix doesn't address?

brandonshensley commented 3 years ago

There is undoubtedly a nicer way to share this, but I've put together a little proof of concept here: https://www.dropbox.com/s/zl64y70bsouco99/polfractens.ipynb?dl=0

Given that this looks reasonable, the current thought is that there is a bug in the code joining the small and large scales in alm-space. Stay tuned.

giuspugl commented 3 years ago

Summary of the analysis performed so far

The tests performed with @brandonshensley are converging to the fact that using a map with variable resolution (e.g. Planck 2015 I map) can be the reason why we see a knee in the power spectra of real quantities as stated above by @seclark. To further simplify the test we therefore proceeded by considering the Planck 2018 unires map which has the downside of being coarser in terms of resolution , 80' , and of being less optimized in the removal of the CIB but as a pro it has the same resolution in each region of the sky. We repeated the same procedure as outlined above and we restricted in particular to the case with P=0, i.e. i = log( I ) The procedure so far i used :

i = log(I)
run Namaster on 80% sky on the i map
fit from the Namaster spectrum the power law to get the small scales ss map
1. get the large scale map with ls= Alm2map(almxfl(i,sigmoid))
2. I’ = exp(ls +ss) we combine the maps ls and ss in step 5 but we can't combine ss and ls together as the two maps are estimated on different domains, the former estimated by fitting the spectra estimated with namaster and the other estimated from pseudo -alm . We fix this issue by running 2., 3. and 4. on fullsky quantities (i.e. using anafast) .

I map at 80'

We transform to the i map and estimate the power spectrum with anafast

With a lower reso, we have to restrict our fitting range in ell=30-100 .The spectral index of the fitted power law is -1.58. A comparison on the map level of the maps with and without small scales injected by extrapolating the power law is shown below :

and the TT power spectrum estimated from the I' map (solid orange) , compared with the input I with 80' reso (solid blue). As you can see the spectra behave as expected. Notice that for the real quantities we can keep estimating the spectra with Namaster on the 80% of the sky.

Modulating small scales with a template map.

We now want to inject a modulation of the small scales similarly to what was done in the PySM2, since so far the amplitude is uniform across the sky.
We therefore : a. smoothed the same I map used above to 5 deg, b. take the log10, c. saturate to 2 all the pixels whose log10 value is >2. d. Minmax rescale the map from 0.1 to 1 . In this way we shrink the dynamic range where we expect the scales to be modulated by the dust intensity profile, A . Below the map produced following the steps from a-d .

Then we follow the steps from 1-5 and just changed the definition in 5. as : I’ = exp(ls +( A*ss )) below you can see a direct comparison of the maps produced with the modulation defined above for a region at ~intermediate Gal. latitudes (same area as the one shown above) Notice how the small scales now are dimmer in the vicinity of low intensity regions and viceversa.

Finally, we want to check whether this modulation affects the power spectrum.

the anafast power spectrum obtained from the logarithmic map (shown in solid orange) do not seem to be affected by the modulation ( compared with the input i map spectrum shown in solid blue) .

As well as the TT power spectrum estimated from the I' map (solid orange) , compared with the input I with 80' reso (solid blue).

Remarks

We propose to use the Planck 2018 GNILC maps with unires both for IQU maps, given the fact that we don't see the inconsistency observed in TT spectra when varres maps are chosen.
We found a modulation technique that seems to be working both at the map and at the power spectrum level ( we can further tweak some of the parameters, e.g. resolution of the template map, saturation scale, etc...)
didn't have time to clean up the notebook will do asap (apologies!)

seclark commented 3 years ago

This sounds really promising @giuspugl! Glad to hear that the tests are converging. From a conversation with @brandonshensley, I am curious now what the histogram of the polarization fraction looks like after the addition of small scales

giuspugl commented 3 years ago

Finalizing the injection of small scales on GNILC dust maps

With the last commit , we show reasonable results in injecting small scales in the poltens formalism on the dust emission.

We use the 2015 GNILC intensity map from the 2nd Planck release, as it encodes less contamination from CIB at 21.8' resolution
for Q and U we adopt maps from the 3rd Planck release as they were optimized for polarization studies at 80' reso.

in this way, we make our modeling of smaller angular scale to fully benefit of the currently available data sets, where Intensity maps have been assessed at good high SNR at smaller resolutions than polarization ones. We summarize the preprocessing steps as it follows:

Monopole subtraction: we subtract the value reported in Section 2.2 of Planck 2018 XII to the I map for the pixels outside the Galactic 080 mask .
We fit the power-law in different multipole ranges as the input maps are at different resolutions
as the injected small scales are at different multipoles for intensity and polarization, we consider 2 different modulation maps
estimate the NaMaster spectra outside 080 mask of the maps before and after the small scale injection
estimate the polarization fraction distribution before and after the small scale injection
eventually re-include the monopole term to the intensity map ( not done yet ...) The updated jupyter notebook is available here

zonca commented 3 years ago

thanks @giuspugl

do you think this is completed? should I start turning this into production code to be merged into PySM3? Or better first produce IQU maps to be shared with the Panex group for further analysis?

giuspugl commented 3 years ago

We decided that we 're ready to port it into PySM3, please let me know if I can be of any help on this!

zonca commented 3 years ago

very good! I'll start to work on this and keep you posted.

zonca commented 3 years ago

@giuspugl should I use the gnilc_gaussian_scales.ipynb notebook inside this pull request or expansion_logpoltens.ipynb in the Gist?

giuspugl commented 3 years ago

@zonca i recommend to use the notebook in the gist here

zonca commented 3 years ago

@giuspugl can you please make /global/cscratch1/sd/giuspugl/workstation/FG_extension/extending_gnilc_dust/ readable so I can compare results?

giuspugl commented 3 years ago

Sure, i made a new directory /global/cscratch1/sd/giuspugl/workstation/FG_extension/extending_gnilc_dust you 'd be able to access .

zonca commented 3 years ago

continues in #82

galsci / pysm