[frozen] Gauss Legendre pixelization implementation of d11

zonca commented 2 years ago

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

zonca commented 2 years ago

Comparison of dx11 and dx11GL

I am debugging a difference between the HEALPix and the GL implementation of dx11. I boiled it down to the logpoltens transform.

If I have:

Same Alm with large and small scales
I transform to pixel space (GL or HEALPix)
Apply log pol tens anti-transform to get back to IQU
Turn GL back to Alm and then to HEALPix
Compare the 2 HEALPix maps, one that went through GL the other one direct

If I skip the logpoltens anti-transform, comparison is fine:

Applying the anti-transform I get:

Gnomview of galactic center:

Full emission is:

Notebook with test: https://gist.github.com/325e796531f8d17acce1d258b0c4752e

Logpoltens to IQU transform

https://github.com/galsci/pysm/blob/f760ad26363dcf4f4ee235a791342699b038d9d9/pysm3/utils/logpoltens.py#L19-L29

@mreineck maybe I am missing something obvious? @giuspugl @seclark @brandonshensley any suggestion?

mreineck commented 2 years ago

Well, logpoltens is a nonlinear transform, so it destroys the band limitedness property of the map - and depending on the pixelization this has different effects on map analysis.

Is it absolutely necessary for the workflow to have a_lm of logpoltens quantities? As long as you stick to transforming IQU maps, everything should be fine. But running logpoltens quantities through map2alm (no matter in which pixelization) feels ill-defined to me.

zonca commented 2 years ago

thanks @mreineck! is it also possible that GL pixels having different areas also impacts the result? Should I weight them before applying the logpoltens anti-transform?

Well, logpoltens is a nonlinear transform, so it destroys the band limitedness property of the map - and depending on the pixelization this has different effects on map analysis.

Is it absolutely necessary for the workflow to have a_lm of logpoltens quantities? As long as you stick to transforming IQU maps, everything should be fine. But running logpoltens quantities through map2alm (no matter in which pixelization) feels ill-defined to me.

We run map2alm on logpoltens quantities but only keep the largest scales. The main issue here, I think, is that we generate the small scales in logpoltens A_lm, so they are band-limited in logpoltens, but when we anti-transform to IQU maps it is now the IQU maps which are not band-limited.

mreineck commented 2 years ago

is it also possible that GL pixels having different areas also impacts the result? Should I weight them before applying the logpoltens anti-transform?

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.

We run map2alm on logpoltens quantities but only keep the largest scales. The main issue here, I think, is that we generate the small scales in logpoltens A_lm, so they are band-limited in logpoltens, but when we anti-transform to IQU maps it is now the IQU maps which are not band-limited.

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.

zonca commented 2 years ago

Spectra of the residual maps

We are comparing:

HEALPix = lopoltens2map(alm2map (small + large scales) )
HEALPix(Alm) = alm2map (map2alm(HEALPix) )
GL= alm2map(GL2alm(lopoltens2map(alm2GL(small + large scales))))

Notebook with the comparison: https://gist.github.com/1b74176a66321aa451efc9b77cc0888a

Clip small scales at 2.5 nside

zonca commented 2 years ago

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:

    rtol = 1e-4
    assert_quantity_allclose(output_healpix[0], output_gl_to_healpix[0], rtol=rtol, atol=5*u.uK_RJ)
    assert_quantity_allclose(output_healpix[1:], output_gl_to_healpix[1:], rtol=rtol, atol=1*u.uK_RJ)

Suggestions?

zonca commented 2 years ago

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,
    )

zonca commented 2 years ago

next let's compare making the smoothing on the templates (as proposed in dx11) or after generating the templates (as is currently done in dx10):

zonca commented 2 years ago

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. It is actually pretty fast to run small scales on the fly. So I suggest we could generate small scales on the fly also for dx10, just fix the seeds and run with 16k Lmax to that we always get the same results. So basically have both dx10 and dx9 that call the machinery of dx11 with fixed seeds and with fixed spectral index for dx9 instead of running the same class of dx1. Then we can provide also dx9gl dx10gl and dx11gl which use Gauss-Legendre pixelization, which is faster and more accurate, but requires ducc0. Let's discuss at next suitable Panexp call.

zonca commented 2 years ago

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. Open to suggestions.

zonca commented 2 years ago

Had a discussion today with @seclark and @delabru, the agreement between GL and HEALPix implementation of dx11 is not good enough, it is just going to be confusing if we publish both types of models.

So the plan is to keep this pull request open and unmerged, then revisit this later on. Anyway, the implementation is complete, the unit tests run and pass (comparing the smoothed maps).

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.

galsci / pysm