-
Notifications
You must be signed in to change notification settings - Fork 23
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
[frozen] Gauss Legendre pixelization implementation of d11 #116
base: main
Are you sure you want to change the base?
Conversation
zonca
commented
Apr 20, 2022
- feat: add ModifiedBlackBodyRealizationGL as a copy of ModifiedBlackBodyRealization
- feat: map2alm and alm2map for Gauss Legendre with test
- utils
- feat: first running implementation of dx11-style realization in GL pixelization
Comparison of dx11 and dx11GLI am debugging a difference between the HEALPix and the GL implementation of If I have:
If I skip the logpoltens anti-transform, comparison is fine: Applying the anti-transform I get: Full emission is: Notebook with test: https://gist.github.com/325e796531f8d17acce1d258b0c4752e Logpoltens to IQU transformpysm/pysm3/utils/logpoltens.py Lines 19 to 29 in f760ad2
@mreineck maybe I am missing something obvious? |
Well, Is it absolutely necessary for the workflow to have a_lm of |
thanks @mreineck! is it also possible that GL pixels having different areas also impacts the result? Should I weight them before applying the
We run |
All the weighting is taken care of by the SHTs internally (otherwise you wouldn't see the really tiny errors when you omit the transform). Other than in Healpix, the GL pixelisation has exactly one set of weights that should be used, and there isn't really any choice.
I see. I don't really have a recommendation here. Even if things look better when using Healpix, I fear that the Healpix pixelization just hides the problems better, but they still exist. |
Spectra of the residual mapsWe are comparing:
Notebook with the comparison: https://gist.github.com/1b74176a66321aa451efc9b77cc0888a Clip small scales at 2.5 nside |
ok, I investigated a bit more. I realized that we always want to smooth maps after we have generated them with the beam of the instrument, so I compared the smoothed maps instead of the modelling maps. Of course, now the difference between ell = 2 nside and 3 nside is not as bad. See notebook here: https://gist.github.com/75a6929e14b95b8f3f729e1a41d68941 Here are the 2 spectra, nside 512 smoothed with 14 arcmin beam. Still the ratio shows residual difference between GL and Healpix: However, doesn't impact much the output maps: So unit tests pass with this tolerance:
Suggestions? |
I might have found a better solution: move the beam smoothing upfront, users can still choose to run with no beam, but the preferred way of running a model would be to provide a beam to the Sky object, then dx11 and dx11gl can directly apply the beam to the small scales, so that the logpoltens map is band limited, now even after the non linear logpoltens transform, there is less power at small scales and the spherical harmonics transforms are more well behaved. Now consistency between dx11 and dx11gl is better: rtol = 1e-4
atol = {353: 20 * u.uK_RJ, 545: 60 * u.uK_RJ, 857: 150 * u.uK_RJ}
assert_quantity_allclose(
output_healpix[0], output_gl_to_healpix[0], rtol=rtol, atol=atol[freq.value]
)
assert_quantity_allclose(
output_healpix[1:],
output_gl_to_healpix[1:],
rtol=rtol,
atol=atol[freq.value] / 10,
) |
what I am concerned about now is dx9 and dx10. They have been generated with significant power at high ell, so the spherical harmonics transforms are introducing artifacts. |
alternatively rerun dx9 and dx10 with presmoothed templates at every pixelization. So for example if we are running at nside 512, we smooth the templates at 2 pixels per beam ~ 14 armin I think. Then people will have to apply the extra smoothing to get to their desired resolution. |
Had a discussion today with @seclark and @delabru, the agreement between GL and HEALPix implementation of So the plan is to keep this pull request open and unmerged, then revisit this later on. However, we still want to better understand the impact of the logpoltens transform on the map, it is useful to have a caveat about this in the paper and in the documentation. It is possible that the logpoltens anti-transform makes the map not band limited, so there is residual power at high ell that could impact analysis. Will open a separate issue to track this. |