#144 Different alms with anafast bwteen Healpix (fortran) and healpy

[Anonymous]

Dear all,

I am Yoshinori.
I got different alms (cls) from the same map with anafast of Healpix build on Fortran and healpy.
Comparison of those cls are attached.
I am being confused because those should be the same in my understanding.
Maps from those alms with alm2map function indicate that the alms from healpy is more plausible than Healpix.
Configurations that I am using is written below. I really appreciate if I can get some reply.

Map: Planck 100GHz map (HFI_SkyMap_100-BPassCorrected_2048_R4.00_full.fits)
Healpix: 3.83
healpy: 1.18.0

Function for healpix
simul_type= 2
nlmax= 6143
theta_cut_deg= 0
iter_order= 3
infile= //HFI_SkyMap_100-BPassCorrected_2048_R4.00_full.fits
outfile= //cl.fits
regression= 0
outfile_alms = /*/alms1.fits

Function for healpy
m1 = m1, h = hp.read_map('/HFI_SkyMap_100-BPassCorrected_2048_R4.00_full.fits', field = (0,1,2,),h=True)
nside = 2048
lmax = nside3-1
mmax = lmax
clp1, almp1 = hp.anafast(m1, lmax=lmax, mmax=mmax, alm=True)

[Martin Reinecke]

Hi Yoshinori,

thank you for submitting this!
From the attached plot it almost looks like the Fortran code has been run without iterating (deviation in TT at low l, noise at high l), even though you specified iter_order=3. Just as an experiment, could you re-run healpy with 0 iterations and compare that to the Healpix results?
Another strange thing is that the TT spectra at low l match very well for even l, but badly for odd ones. I'm pretty sure this is a hint at the source of the problem, but I really have no concrete idea at the moment ... I'll think some more.

Best regards,
Martin

[Anonymous]

Dear Martin,

Thank you for you reply.
I attached plots you requests: "Healpix_iter3 vs Healpy_ite0" and "Healpix_iter3 vs Healpix_ite0" .
I confirmed Fortran code output from iter3 and iter0 are different, which indicate it worked somehow.
On the other hand, healpy output with iter0 was not different than healpy output with iter3.
Healpix_iter3 vs Healpy_ite0 has still difference.
I also attached map comparizon between INPUT map and OUTPUT map calculated with healpy.synfast by using the output alms in Fortran code.
I hope this would be helpful.

For more reference,
I confirmed healpy anafast vs healpix anafast using simulated maps using cl.fits that stored in healpix package (in test directory). Configuration is slightly difference; I used Nside = 128. I found no significant difference like what I originally shown although there is slightly difference in terms of alms.

I am confused both anafast (healpy and Fortran code) works well at some maps but do not work at another maps.

If there is things what you want me to do, I am happy to do that. Please let me know.
Thank you for your consideration.

Best Regards,
Yoshinori

[Anonymous]

Do you see my post above?

[Eric Hivon]

Hi Yoshinori, I see the difference maps, but not the the plot of C(l) with various iter values.
The Planck map you are analyzing is real data, and not necessarily band limited (ie, it has power at very small scales, at least because of noise), whereas the simulated map you read from the Healpix package is band limited by construction, because its lmax was set to 2xNside.
Beware that healpy and Healpix may not treat the monopole and dipole the same way.
I think that by default healpy with regress out the best fit dipole and monopole before computing the alm, while you have to request for it in Healpix.

[Anonymous]

Thank you for your reply.
It seems I can attach only one figure.
This is Healpix_iter3 vs Healpy_ite0.

[Anonymous]

This is "Healpix_iter3 vs Healpy_ite0"

[Anonymous]

The Planck map you are analyzing is real data, and not necessarily band limited (ie, it has power at very small scales, at least because of noise), whereas the simulated map you read from the Healpix package is band limited by construction, because its lmax was set to 2xNside.

Yeah, I read such kind of documents and understood it, but thank you for reminding.
We should get the same alms when a input map is the same. That's why I am confused.

Beware that healpy and Healpix may not treat the monopole and dipole the same way.
I think that by default healpy with regress out the best fit dipole and monopole before computing the alm, while you have to request for it in Healpix.

But this looks the reason. So does it mean Cpp codes and Fortran codes are different for monopole diople treatment? Since I could not the documents, I would happy to let me know where it is.
I checked Fortran code with regression = 2 (theta_cut_deg = 0deg), which remove monopole and dipole, but there is still difference as shown in attached file although the cls gots much closer.
Is there any way to check the difference of monopole/dipole treatment between healpy and healpix?

Best Regards,
Yoshinori

[Anonymous]

Any other comments or suggestions for the reason or debug method that is worth doing?

Thanks,
Yoshinori

[Eric Hivon]

HI Yoshinori, thanks for the plots.
Just to be sure, those latest plots have Healpix run with the dipole+monopole regression AND iter=3 ?
The regression only affects the T map, and therefore only the TT spectrum.
The discrepancy seen in EE and BB is not related to that.
And did you try running the fortran code in double precision ?