diff --git a/docs/modules/healpix.rst b/docs/modules/healpix.rst index a551a0775..1a38bd79e 100644 --- a/docs/modules/healpix.rst +++ b/docs/modules/healpix.rst @@ -16,7 +16,6 @@ This module should be imported manually if used outside of GLASS:: --------------------- .. autofunction:: ang2pix -.. autofunction:: ang2vec .. autofunction:: npix2nside .. autofunction:: nside2npix .. autofunction:: randang @@ -24,11 +23,6 @@ This module should be imported manually if used outside of GLASS:: ``healpy`` Functions -------------------- -.. autofunction:: alm2map -.. autofunction:: alm2map_spin -.. autofunction:: almxfl -.. autofunction:: get_nside -.. autofunction:: map2alm .. autofunction:: pixwin .. autofunction:: query_strip .. automethod:: Rotator.rotate_map_pixel diff --git a/examples/2-advanced/legacy-mode.ipynb b/examples/2-advanced/legacy-mode.ipynb index 5743dc00f..59b9cabc7 100644 --- a/examples/2-advanced/legacy-mode.ipynb +++ b/examples/2-advanced/legacy-mode.ipynb @@ -494,7 +494,7 @@ { "data": { "text/plain": [ - "np.False_" + "False" ] }, "execution_count": null, @@ -537,7 +537,7 @@ { "data": { "text/plain": [ - "np.True_" + "True" ] }, "execution_count": null, diff --git a/glass/fields.py b/glass/fields.py index 8e59c217f..c15ec5151 100644 --- a/glass/fields.py +++ b/glass/fields.py @@ -394,13 +394,7 @@ def _generate_grf( alm = xpx.at(alm)[:n].set(xp.real(alm[:n]) + xp.imag(alm[:n]) + 0j) # ty: ignore[not-subscriptable] # transform alm to maps - yield glass.harmonics.inverse_transform( - alm, - lmax=n - 1, - nside=nside, - pixwin=False, - pol=False, - ) + yield glass.harmonics.inverse_transform(alm, lmax=n - 1, nside=nside) def getcl( diff --git a/glass/harmonics.py b/glass/harmonics.py index 05f6ec951..09a81d503 100644 --- a/glass/harmonics.py +++ b/glass/harmonics.py @@ -11,7 +11,7 @@ import array_api_compat -import glass.healpix as hp +from glass.healpix import _alm2map, _alm2map_spin, _map2alm if TYPE_CHECKING: from glass._types import ComplexArray, FloatArray @@ -60,8 +60,6 @@ def inverse_transform( *, lmax: int, nside: int, - pixwin: bool = False, - pol: bool = True, spin: typing.Literal[0] = 0, ) -> FloatArray: # returns a real scalar map @@ -74,8 +72,6 @@ def inverse_transform( *, lmax: int, nside: int, - pixwin: bool = False, - pol: bool = True, spin: typing.Literal[1, 2], ) -> ComplexArray: # returns a complex map for spin 1 or 2 @@ -88,21 +84,17 @@ def inverse_transform( *, lmax: int, nside: int, - pixwin: bool = False, - pol: bool = True, spin: int, ) -> FloatArray | ComplexArray: # returns a real or complex map depending on the spin ... -def inverse_transform( # noqa: PLR0913 +def inverse_transform( alm: ComplexArray, *, lmax: int, nside: int, - pixwin: bool = False, - pol: bool = True, spin: int = 0, ) -> FloatArray | ComplexArray: """ @@ -120,11 +112,6 @@ def inverse_transform( # noqa: PLR0913 The maximum multipole of the spherical harmonic transform. nside The nside parameter of the output map. - pixwin - Whether to apply the pixel window function. Has no effect for non-zero - spin. - pol - Whether to compute polarization. Has no effect for non-zero spin. spin Spin of the output map. Zero produces a real scalar map; non-zero spin produces a complex map whose real and imaginary parts are the two @@ -138,9 +125,9 @@ def inverse_transform( # noqa: PLR0913 xp = alm.__array_namespace__() if spin == 0: - return hp.alm2map(alm, nside, lmax=lmax, pixwin=pixwin, pol=pol) + return _alm2map(alm, nside, lmax=lmax) - maps = hp.alm2map_spin([alm, xp.zeros_like(alm)], nside, spin, lmax) + maps = _alm2map_spin([alm, xp.zeros_like(alm)], nside, spin, lmax) return maps[0] + 1j * maps[1] @@ -148,8 +135,6 @@ def transform( maps: FloatArray, *, lmax: int, - pol: bool = True, - use_pixel_weights: bool = False, ) -> ComplexArray: """ Compute the spherical harmonic transform of a map. @@ -160,14 +145,11 @@ def transform( The real-space map to transform. lmax The maximum multipole of the spherical harmonic transform. - pol - Whether to compute polarization. - use_pixel_weights - Whether to use pixel weights in the transform. Returns ------- The spherical harmonic coefficients resulting from the transform. """ - return hp.map2alm(maps, lmax=lmax, pol=pol, use_pixel_weights=use_pixel_weights) + # the only time GLASS used _map2alm was with pixel weights + return _map2alm(maps, lmax=lmax, use_pixel_weights=True) diff --git a/glass/healpix.py b/glass/healpix.py index 90597e92c..bb6756f83 100644 --- a/glass/healpix.py +++ b/glass/healpix.py @@ -38,13 +38,11 @@ def _get_healpy_datapath() -> str | None: @numpy_fallback -def alm2map( +def _alm2map( alms: ComplexArray | Sequence[ComplexArray], nside: int, *, lmax: int | None = None, - pixwin: bool = False, - pol: bool = True, ) -> FloatArray: """ Computes a HEALPix map given the alm. @@ -57,10 +55,6 @@ def alm2map( The nside of the output map. lmax Explicitly define lmax. - pixwin - Smooth the alm using the pixel window functions. - pol - If True, assumes input alms are TEB. Returns ------- @@ -68,17 +62,11 @@ def alm2map( containing the T, Q, and U maps for polarized input. """ - return healpy.alm2map( - alms, - nside, - lmax=lmax, - pixwin=pixwin, - pol=pol, - ) + return healpy.alm2map(alms, nside, lmax=lmax) @numpy_fallback -def alm2map_spin( +def _alm2map_spin( alms: Sequence[FloatArray], nside: int, spin: int, @@ -107,30 +95,6 @@ def alm2map_spin( return list(outputs) -@numpy_fallback -def almxfl( - alm: ComplexArray, - fl: FloatArray, -) -> ComplexArray: - """ - Multiply alm by a function of l. The function is assumed to be zero where - not defined. - - Parameters - ---------- - alm - The alm to multiply. - fl - The function (at l=0..fl.shape[0]-1) by which alm must be multiplied. - - Returns - ------- - The new alm array with the function of l applied. - - """ - return healpy.almxfl(alm, fl) - - @numpy_fallback def ang2pix( nside: int, @@ -167,61 +131,10 @@ def ang2pix( @numpy_fallback -def ang2vec( - theta: float | FloatArray, - phi: float | FloatArray, - *, - lonlat: bool = False, -) -> tuple[FloatArray, FloatArray, FloatArray]: - """ - Convert angles to 3D position vector. - - Parameters - ---------- - theta - Angular coordinates of a point on the sphere. - phi - Angular coordinates of a point on the sphere. - lonlat - If True, automatically adjust latitudes to be within [-90, 90] range. - - Returns - ------- - A normalised 3-vector pointing in the same direction as ``ang``. - - """ - x, y, z = healpix.ang2vec( - theta, - phi, - lonlat=lonlat, - ) - return x, y, z - - -@numpy_fallback -def get_nside(m: FloatArray) -> int: - """ - Return the nside of the given map. - - Parameters - ---------- - m - The map to get the nside from. - - Returns - ------- - The HEALPix nside parameter of the map. - - """ - return int(healpy.get_nside(m)) - - -@numpy_fallback -def map2alm( +def _map2alm( maps: FloatArray | Sequence[FloatArray], *, lmax: int | None = None, - pol: bool = True, use_pixel_weights: bool = False, ) -> ComplexArray: """ @@ -243,9 +156,6 @@ def map2alm( The input map or a list of n input maps. Must be in ring ordering. lmax Maximum l of the power spectrum. - pol - If True, assumes input maps are TQU. If False, only the temperature map - is considered. use_pixel_weights If True, use pixel by pixel weighting, healpy will automatically download the weights, if needed. @@ -259,7 +169,6 @@ def map2alm( maps, datapath=_get_healpy_datapath(), lmax=lmax, - pol=pol, use_pixel_weights=use_pixel_weights, ) diff --git a/glass/lensing.py b/glass/lensing.py index 4108269c4..6a7cf2d6e 100644 --- a/glass/lensing.py +++ b/glass/lensing.py @@ -289,17 +289,12 @@ def from_convergence( # noqa: PLR0913 return () # get the NSIDE parameter - nside = hp.get_nside(kappa) + nside = hp.npix2nside(kappa.shape[-1]) if lmax is None: lmax = 3 * nside - 1 # compute alm - alm = glass.harmonics.transform( - kappa, - lmax=lmax, - pol=False, - use_pixel_weights=True, - ) + alm = glass.harmonics.transform(kappa, lmax=lmax) # mode number; all conversions are factors of this. # Must be float to allow division later @@ -315,7 +310,7 @@ def from_convergence( # noqa: PLR0913 lambda a: xp.divide(-2.0, a * (a + 1)), fill_value=0, ) - alm = hp.almxfl(alm, fl) + alm = glass.harmonics.multalm(alm, fl) # if potential is requested, compute map and add to output if potential: @@ -329,7 +324,7 @@ def from_convergence( # noqa: PLR0913 # compute deflection alms in place fl = xp.sqrt(ell * (ell + 1)) # missing spin-1 pixel window function here - alm = hp.almxfl(alm, fl) + alm = glass.harmonics.multalm(alm, fl) # if deflection is requested, compute spin-1 maps and add to output if deflection: @@ -352,7 +347,7 @@ def from_convergence( # noqa: PLR0913 if discretized: pw0, pw2 = hp.pixwin(nside, lmax=lmax, pol=True, xp=xp) fl *= pw2 / pw0 - alm = hp.almxfl(alm, fl) + alm = glass.harmonics.multalm(alm, fl) # transform to shear maps gamma = glass.harmonics.inverse_transform(alm, lmax=lmax, nside=nside, spin=2) diff --git a/tests/core/test_harmonics.py b/tests/core/test_harmonics.py index c5181a62d..db15ad3aa 100644 --- a/tests/core/test_harmonics.py +++ b/tests/core/test_harmonics.py @@ -58,7 +58,7 @@ def test_inverse_transform_healpy(xp: ModuleType) -> None: nside = 1 lmax = 2 - expected = hp.alm2map(alm, nside=nside, lmax=lmax) + expected = hp._alm2map(alm, nside=nside, lmax=lmax) actual = glass.harmonics.inverse_transform(alm, lmax=lmax, nside=nside) xpx.testing.assert_equal(actual, expected) @@ -72,7 +72,7 @@ def test_inverse_transform_healpy_spin( lmax = 2 nside = 1 - result = hp.alm2map_spin([alm, xp.zeros_like(alm)], nside, spin, lmax) + result = hp._alm2map_spin([alm, xp.zeros_like(alm)], nside, spin, lmax) expected = result[0] + 1j * result[1] actual = glass.harmonics.inverse_transform(alm, lmax=lmax, nside=nside, spin=spin) xpx.testing.assert_equal(actual, expected) @@ -82,6 +82,11 @@ def test_transform_healpy(xp: ModuleType) -> None: maps = xp.asarray([1.0] * 12) lmax = 2 - expected = hp.map2alm(maps, lmax=lmax) + expected = hp._map2alm( + maps, + lmax=lmax, + # using pixel weights is the default in the transform + use_pixel_weights=True, + ) actual = glass.harmonics.transform(maps, lmax=lmax) xpx.testing.assert_equal(actual, expected) diff --git a/tests/core/test_healpix.py b/tests/core/test_healpix.py index 3073014f9..6a371c93c 100644 --- a/tests/core/test_healpix.py +++ b/tests/core/test_healpix.py @@ -19,60 +19,34 @@ from tests.fixtures.helper_classes import HealpixInputs -@pytest.mark.parametrize( - ("pixwin", "pol"), - [ - (False, False), - (False, True), - (True, False), - (True, True), - ], -) -def test_alm2map_individual( # noqa: PLR0913, PLR0917 +def test_alm2map_individual( healpix_inputs: type[HealpixInputs], np: ModuleType, - pixwin: bool, # noqa: FBT001 - pol: bool, # noqa: FBT001 urng: UnifiedGenerator, xp: ModuleType, ) -> None: - """Compare ``glass.healpix.alm2map`` against ``healpy.alm2map``.""" + """Compare ``glass.healpix._alm2map`` against ``healpy.alm2map``.""" alm = healpix_inputs.alm(rng=urng) old = healpy.alm2map( np.asarray(alm), healpix_inputs.nside, lmax=healpix_inputs.lmax, - pixwin=pixwin, - pol=pol, ) - new = hp.alm2map( + new = hp._alm2map( alm, healpix_inputs.nside, lmax=healpix_inputs.lmax, - pixwin=pixwin, - pol=pol, ) xpx.testing.assert_equal(xp.asarray(old), new) -@pytest.mark.parametrize( - ("pixwin", "pol"), - [ - (False, False), - (False, True), - (True, False), - (True, True), - ], -) -def test_alm2map_sequence( # noqa: PLR0913, PLR0917 +def test_alm2map_sequence( healpix_inputs: type[HealpixInputs], np: ModuleType, - pixwin: bool, # noqa: FBT001 - pol: bool, # noqa: FBT001 urng: UnifiedGenerator, xp: ModuleType, ) -> None: - """Compare ``glass.healpix.alm2map`` against ``healpy.alm2map``.""" + """Compare ``glass.healpix._alm2map`` against ``healpy.alm2map``.""" alm = healpix_inputs.alm(rng=urng) blm = healpix_inputs.alm(rng=urng) clm = healpix_inputs.alm(rng=urng) @@ -80,15 +54,11 @@ def test_alm2map_sequence( # noqa: PLR0913, PLR0917 [np.asarray(alm), np.asarray(blm), np.asarray(clm)], healpix_inputs.nside, lmax=healpix_inputs.lmax, - pixwin=pixwin, - pol=pol, ) - new = hp.alm2map( + new = hp._alm2map( [alm, blm, clm], healpix_inputs.nside, lmax=healpix_inputs.lmax, - pixwin=pixwin, - pol=pol, ) xpx.testing.assert_equal(xp.asarray(old), new) @@ -100,7 +70,7 @@ def test_alm2map_spin( urng: UnifiedGenerator, xp: ModuleType, ) -> None: - """Compare ``glass.healpix.alm2map_spin`` against ``healpy.alm2map_spin``.""" + """Compare ``glass.healpix._alm2map_spin`` against ``healpy.alm2map_spin``.""" alm = healpix_inputs.alm(rng=urng) blm = healpix_inputs.alm(rng=urng) old = healpy.alm2map_spin( @@ -109,26 +79,13 @@ def test_alm2map_spin( spin, healpix_inputs.lmax, ) - new = hp.alm2map_spin([alm, blm], healpix_inputs.nside, spin, healpix_inputs.lmax) + new = hp._alm2map_spin([alm, blm], healpix_inputs.nside, spin, healpix_inputs.lmax) assert type(old) is type(new) assert len(old) == len(new) for i in range(len(old)): xpx.testing.assert_equal(xp.asarray(old[i]), new[i]) -def test_almxfl( - healpix_inputs: type[HealpixInputs], - urng: UnifiedGenerator, - xp: ModuleType, -) -> None: - """Compare ``glass.healpix.almxfl`` against ``healpy.almxfl``.""" - alm = healpix_inputs.alm(rng=urng) - fl = healpix_inputs.fl(rng=urng) - old = healpy.almxfl(alm, fl) - new = hp.almxfl(alm, fl) - xpx.testing.assert_equal(xp.asarray(old), new) - - @pytest.mark.parametrize( ("lonlat", "max_phi", "max_theta"), [ @@ -152,139 +109,84 @@ def test_ang2pix( # noqa: PLR0913,PLR0917 xpx.testing.assert_equal(xp.asarray(old), new) -@pytest.mark.parametrize( - ("lonlat", "max_phi", "max_theta"), - [ - (False, math.pi, math.pi / 2), - (True, 180, 90), - ], -) -def test_ang2vec( # noqa: PLR0913,PLR0917 - healpix_inputs: type[HealpixInputs], - lonlat: bool, # noqa: FBT001 - max_phi: float, - max_theta: float, - urng: UnifiedGenerator, - xp: ModuleType, -) -> None: - """Compare ``glass.healpix.ang2vec`` against ``healpix.ang2vec``.""" - thetas = healpix_inputs.longitudes(max_theta, rng=urng) - phis = healpix_inputs.latitudes(max_phi, rng=urng) - old = healpix.ang2vec(thetas, phis, lonlat=lonlat) - new = hp.ang2vec(thetas, phis, lonlat=lonlat) - assert type(old) is type(new) - assert len(old) == len(new) - for i in range(len(old)): - xpx.testing.assert_equal(xp.asarray(old[i]), new[i]) - - -def test_get_nside( - healpix_inputs: type[HealpixInputs], - np: ModuleType, - urng: UnifiedGenerator, -) -> None: - """Compare ``glass.healpix.get_nside`` against ``healpy.get_nside``.""" - kappa = healpix_inputs.kappa(rng=urng) - assert healpy.get_nside(np.asarray(kappa)) == hp.get_nside(kappa) - - @pytest.mark.usefixtures("_add_healpy_datapath_to_env") def test_map2alm_with_pulled_data( healpix_inputs: type[HealpixInputs], urng: UnifiedGenerator, ) -> None: - """Tests running map2alm offline works as expected.""" + """Tests running _map2alm offline works as expected.""" npix = healpy.nside2npix(32) kappa = urng.normal(size=npix) - result = hp.map2alm( + result = hp._map2alm( kappa, lmax=healpix_inputs.lmax, - pol=False, use_pixel_weights=True, ) assert result.shape == (78,) -@pytest.mark.parametrize("pol", [False, True]) def test_map2alm_with_pulled_data_wrong_path( invalid_healpy_datapath: str, healpix_inputs: type[HealpixInputs], - pol: bool, # noqa: FBT001 urng: UnifiedGenerator, ) -> None: - """Tests running map2alm offline incorrectly doesn't fallback to a HTTP request.""" + """Tests running _map2alm offline incorrectly doesn't fallback to a HTTP request.""" npix = healpy.nside2npix(32) kappa = urng.normal(size=npix) with pytest.raises( ValueError, match=f"Healpy datapath not found at '{invalid_healpy_datapath}", ): - hp.map2alm( + hp._map2alm( kappa, lmax=healpix_inputs.lmax, - pol=pol, use_pixel_weights=True, ) -@pytest.mark.parametrize("pol", [True, False]) @pytest.mark.parametrize("use_pixel_weights", [True, False]) -def test_map2alm_individual( # noqa: PLR0913, PLR0917 +def test_map2alm_individual( healpix_inputs: type[HealpixInputs], np: ModuleType, - pol: bool, # noqa: FBT001 urng: UnifiedGenerator, use_pixel_weights: bool, # noqa: FBT001 xp: ModuleType, ) -> None: - """Compare ``glass.healpix.map2alm`` against ``healpy.map2alm``.""" + """Compare ``glass.healpix._map2alm`` against ``healpy.map2alm``.""" kappa = healpix_inputs.kappa(rng=urng) old = healpy.map2alm( np.asarray(kappa), lmax=healpix_inputs.lmax, - pol=pol, use_pixel_weights=use_pixel_weights, ) - new = hp.map2alm( + new = hp._map2alm( kappa, lmax=healpix_inputs.lmax, - pol=pol, use_pixel_weights=use_pixel_weights, ) xpx.testing.assert_equal(xp.asarray(old), new) -@pytest.mark.parametrize( - ("pol", "use_pixel_weights"), - [ - (False, False), - (False, True), - (True, False), - (True, True), - ], -) -def test_map2alm_sequence( # noqa: PLR0913, PLR0917 +@pytest.mark.parametrize("use_pixel_weights", [True, False]) +def test_map2alm_sequence( healpix_inputs: type[HealpixInputs], np: ModuleType, - pol: bool, # noqa: FBT001 urng: UnifiedGenerator, use_pixel_weights: bool, # noqa: FBT001 xp: ModuleType, ) -> None: - """Compare ``glass.healpix.map2alm`` against ``healpy.map2alm``.""" + """Compare ``glass.healpix._map2alm`` against ``healpy.map2alm``.""" kappa1 = healpix_inputs.kappa(rng=urng) kappa2 = healpix_inputs.kappa(rng=urng) kappa3 = healpix_inputs.kappa(rng=urng) old = healpy.map2alm( [np.asarray(kappa1), np.asarray(kappa2), np.asarray(kappa3)], lmax=healpix_inputs.lmax, - pol=pol, use_pixel_weights=use_pixel_weights, ) - new = hp.map2alm( + new = hp._map2alm( [kappa1, kappa2, kappa3], lmax=healpix_inputs.lmax, - pol=pol, use_pixel_weights=use_pixel_weights, ) xpx.testing.assert_equal(xp.asarray(old), new)