Skip to content
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

Draft
wants to merge 16 commits into
base: main
Choose a base branch
from

Conversation

zonca
Copy link
Member

@zonca 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

@zonca
Copy link
Member Author

zonca commented Apr 25, 2022

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:
image
image
image

Applying the anti-transform I get:

image
image
image

Gnomview of galactic center:
image

Full emission is:

image

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

Logpoltens to IQU transform

def log_pol_tens_to_map(log_pol_tens):
"""Transform from Log pol tensor formalism to IQU
See the notebook about preprocessing GNILC maps in
the documentation of PySM 3 for more details"""
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

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

@mreineck
Copy link

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
Copy link
Member Author

zonca commented Apr 26, 2022

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
Copy link

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
Copy link
Member Author

zonca commented May 2, 2022

Spectra of the residual maps

image

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

image

@zonca
Copy link
Member Author

zonca commented May 19, 2022

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.

image

Still the ratio shows residual difference between GL and Healpix:

image

However, doesn't impact much the output maps:

image

image

image

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
Copy link
Member Author

zonca commented May 19, 2022

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
Copy link
Member Author

zonca commented May 22, 2022

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

image
image

@zonca
Copy link
Member Author

zonca commented May 22, 2022

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
Copy link
Member Author

zonca commented May 22, 2022

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
Copy link
Member Author

zonca commented May 25, 2022

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.

@zonca zonca changed the title Gauss Legendre pixelization implementation of d11 [frozen] Gauss Legendre pixelization implementation of d11 May 25, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

2 participants