simonsobs / pixell

A rectangular pixel map manipulation and harmonic analysis library derived from Sigurd Naess' enlib.
Other
42 stars 32 forks source link

Shape flexibility lost after ducc0 #239

Closed msyriac closed 1 year ago

msyriac commented 1 year ago

Description

curvedsky.alm2map (and possibly other related functions) have lost some flexibility in handling the shapes of inputs. For example, if alm is (nelem,) shaped and the output enmap has shape (1,Ny,Nx), then it returns the following error:

Traceback (most recent call last):
  File "/home/msyriac/repos/falafel/tests/reppixell.py", line 11, in <module>
    cs.alm2map(alm,enmap.empty(shape,wcs))
  File "/home/msyriac/.local/lib/python3.10/site-packages/pixell/curvedsky.py", line 150, in alm2map
    return alm2map_2d(alm, map, ainfo=ainfo, minfo=minfo, spin=spin, deriv=deriv, copy=copy,
  File "/home/msyriac/.local/lib/python3.10/site-packages/pixell/curvedsky.py", line 704, in alm2map_2d
    alm2map_raw_2d(alm[I], tmap, ainfo=ainfo, spin=spin, deriv=deriv, nthread=nthread, verbose=verbose, adjoint=adjoint)
  File "/home/msyriac/.local/lib/python3.10/site-packages/pixell/curvedsky.py", line 837, in alm2map_raw_2d
    alm_full, map_full, ainfo, nthread = prepare_raw(alm, map, ainfo=ainfo, deriv=deriv, nthread=nthread)
  File "/home/msyriac/.local/lib/python3.10/site-packages/pixell/curvedsky.py", line 1365, in prepare_raw
    assert map.shape[:-pixdims] == alm.shape[:-1], "map and alm must agree on pre-dimensions"
AssertionError: map and alm must agree on pre-dimensions

What I Did

The following code will work in v0.17.0 but not in v0.21.0 after the ducc0 transition.

from pixell import enmap,curvedsky as cs,utils as u                                                                                                                          
import numpy as np                                                                                                                                                           
import healpy as hp                                                                                                                                                          

shape,wcs = enmap.fullsky_geometry(res=16.0*u.arcmin,proj='car')                                                                                                             
lmax = 500                                                                                                                                                                   
alm = np.ones(hp.Alm.getsize(lmax))                                                                                                                                          
shape = (1,)+shape                                                                                                                                                           
print(alm.shape)                                                                                                                                                             
print(shape)                                                                                                                                                                 
cs.alm2map(alm,enmap.empty(shape,wcs))
amaurea commented 1 year ago

Fixed in latest version.