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.

Sign up for GitHub

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

Jump to bottom

lensing.lens_map_curved sometimes crashes for single-precision maps #283

Open
zatkins2 opened this issue Jan 9, 2025 · 1 comment · May be fixed by #284
Open

lensing.lens_map_curved sometimes crashes for single-precision maps #283

zatkins2 opened this issue Jan 9, 2025 · 1 comment · May be fixed by #284

Comments

@zatkins2
Copy link
Contributor

zatkins2 commented Jan 9, 2025 •
edited
Loading

  • pixell version: 0.28.0
  • numpy version: 1.26.4
  • Python version: 3.10.16
  • Operating System: rhel 8.10

Description

lensing.lens_map_curved crashes when a single-precision map is requested, given some input alms but not others.

What I Did

Here's an example script that captures the bug.

from pixell import enmap, curvedsky, utils, lensing
import camb
import numpy as np

# make a geometry to put maps onto without aliasing
shape, wcs = enmap.fullsky_geometry(res=4*utils.arcmin, variant='fejer1')
print(shape, wcs)

# get cmb and lensing spectra
pars = camb.set_params(H0=67.5, ombh2=0.022, omch2=0.122, mnu=0.06, omk=0, tau=0.06,
                       As=2e-9, ns=0.965, halofit_version='mead', lmax=4000)
pars.set_for_lmax(4000, lens_potential_accuracy=4)

# calculate results for these parameters
results = camb.get_results(pars)
powers = results.get_cmb_power_spectra(pars, CMB_unit='muK', raw_cl=True)
TT = powers["unlensed_scalar"][:, 0]
PP = powers["lens_potential"][:, 0]
ell = np.arange(len(TT))

# get the lensed cmb map, double prec works
dtype = np.float64

phi_alm = curvedsky.rand_alm(ps=PP, lmax=2000, seed=4)
cmb_alm = curvedsky.rand_alm(ps=TT, lmax=2000, seed=5)
lensedcmb = lensing.lens_map_curved(shape, wcs, phi_alm, cmb_alm, dtype=dtype, verbose=True)[0]
print(lensedcmb.geometry, lensedcmb.dtype)

# get the lensed cmb map, double prec works
phi_alm = curvedsky.rand_alm(ps=PP, lmax=2000, seed=3)
cmb_alm = curvedsky.rand_alm(ps=TT, lmax=2000, seed=4)
lensedcmb = lensing.lens_map_curved(shape, wcs, phi_alm, cmb_alm, dtype=dtype, verbose=True)[0]
print(lensedcmb.geometry, lensedcmb.dtype)

# get the lensed cmb map, single prec works for some maps (note `alms` are in double)
dtype = np.float32

phi_alm = curvedsky.rand_alm(ps=PP, lmax=2000, seed=4)
cmb_alm = curvedsky.rand_alm(ps=TT, lmax=2000, seed=5)
lensedcmb = lensing.lens_map_curved(shape, wcs, phi_alm, cmb_alm, dtype=dtype, verbose=True)[0]
print(lensedcmb.geometry, lensedcmb.dtype)

# get the lensed cmb map, single prec crashes crashes for others (note `alms` are in double)
phi_alm = curvedsky.rand_alm(ps=PP, lmax=2000, seed=3)
cmb_alm = curvedsky.rand_alm(ps=TT, lmax=2000, seed=4)
lensedcmb = lensing.lens_map_curved(shape, wcs, phi_alm, cmb_alm, dtype=dtype, verbose=True)[0]
print(lensedcmb.geometry, lensedcmb.dtype)

with the following output

(2700, 5400) car:{cdelt:[-0.06667,0.06667],crval:[0.03333,0],crpix:[2700.50,1350.50]}
Computing grad map
Computing observed coordinates
Computing alpha map
Computing lensed map
((2700, 5400), car:{cdelt:[-0.06667,0.06667],crval:[0.03333,0],crpix:[2700.50,1350.50]}) float64
Computing grad map
Computing observed coordinates
Computing alpha map
Computing lensed map
((2700, 5400), car:{cdelt:[-0.06667,0.06667],crval:[0.03333,0],crpix:[2700.50,1350.50]}) float64
Computing grad map
Computing observed coordinates
Computing alpha map
Computing lensed map
((2700, 5400), car:{cdelt:[-0.06667,0.06667],crval:[0.03333,0],crpix:[2700.50,1350.50]}) float32
Computing grad map
Computing observed coordinates
Computing alpha map
/home/zatkins/.conda/envs/pixell-test/lib/python3.10/site-packages/pixell/lensing.py:343: RuntimeWarning: invalid value encountered in sqrt
  osint  = (1-ocost**2)**0.5
/home/zatkins/.conda/envs/pixell-test/lib/python3.10/site-packages/pixell/lensing.py:344: RuntimeWarning: invalid value encountered in arcsin
  ophi   = ipos[1] + np.arcsin(sind*grad[1]/osint)
/home/zatkins/.conda/envs/pixell-test/lib/python3.10/site-packages/pixell/lensing.py:346: RuntimeWarning: invalid value encountered in arccos
  return np.array([np.arccos(ocost), ophi]), None
Computing lensed map
Segmentation fault (core dumped)

In the above, if I change the first single-precision case to the following instead

# get the lensed cmb map, single prec works for some maps (note `alms` are now single)
dtype = np.float32

phi_alm = curvedsky.rand_alm(ps=PP, lmax=2000, seed=4, dtype=np.complex64)
cmb_alm = curvedsky.rand_alm(ps=TT, lmax=2000, seed=5, dtype=np.complex64)
lensedcmb = lensing.lens_map_curved(shape, wcs, phi_alm, cmb_alm, dtype=dtype, verbose=True)[0]
print(lensedcmb.geometry, lensedcmb.dtype)

then, instead of working as before, now I get

corrupted double-linked list
Aborted (core dumped)

before even Computing grad map is printed.

Perhaps this is just because of floating-point errors at single-precision? For instance, lenspyx will raise an exception if single-precision lensing is attempted.
@zatkins2
Copy link
Contributor Author

Actually I've observed similar RuntimeWarnings with a crash for double-precision too, but unfortunately the seed was None so I can't reproduce it for you easily

@zatkins2 zatkins2 linked a pull request Jan 14, 2025 that will close this issue
Open
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

Successfully merging a pull request may close this issue.

Lens map curved with lenspyx simonsobs/pixell
1 participant
@zatkins2