Skip to content

SHT beam convolution#114

Open
ArtemBasyrov wants to merge 9 commits intomainfrom
sht_convolution
Open

SHT beam convolution#114
ArtemBasyrov wants to merge 9 commits intomainfrom
sht_convolution

Conversation

@ArtemBasyrov
Copy link
Copy Markdown

Added sht operators for jax_healpy map2alm and alm2map. Added two beam operators, they support convolving all Stokes parameters and all frequencies with one b_ell function (harmonic representation of a symmetric beam); all Stokes parameters with different beams for different frequencies (both are BeamOperator); and each Stokes parameter, each frequency with its own beam in BeamOperatorIQU which takes beam_fl as StokesIQU. All new functions come with tests, and custom reduction rules.

In order to run requires:

  1. @ASKabalan branch of jax_healpy
  2. Version of s2fft with custom_transpose rule implemented, which is required for linear_transpose and transpose() capabilities of furax. It's available on my fork s2fft

Copy link
Copy Markdown

Copilot AI left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pull request overview

Adds spherical-harmonic-transform (SHT) operators and beam-convolution operators to furax, along with algebraic reduction rules and tests to simplify operator compositions.

Changes:

  • Introduce Map2Alm / Alm2Map linear operators backed by jax_healpy, plus an SHTRule reduction.
  • Add BeamOperator / BeamOperatorIQU for symmetric-beam convolution in harmonic space, plus beam composition reduction rules.
  • Add unit tests covering shape/dtype behavior, transpose/inverse expectations, and reduction behavior.

Reviewed changes

Copilot reviewed 6 out of 6 changed files in this pull request and generated 7 comments.

Show a summary per file
File Description
src/furax/math/sht.py New SHT operators and SHTRule for simplifying Map2Alm @ Alm2Map.
src/furax/obs/operators/_beam_operator.py New harmonic-space beam operators and beam composition reduction rules.
src/furax/obs/operators/__init__.py Re-export beam operators from the obs operators package.
src/furax/obs/__init__.py Re-export beam operators from the top-level obs namespace.
tests/math/test_sht.py New tests for SHT operator behavior and SHTRule.
tests/obs/operators/test_beam_operator.py New tests for beam operators and beam reduction rules.

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

Comment thread src/furax/math/sht.py Outdated
Comment thread src/furax/math/sht.py
Comment thread src/furax/math/sht.py Outdated
Comment thread src/furax/math/sht.py Outdated
Comment thread src/furax/obs/operators/_beam_operator.py Outdated
Comment thread src/furax/obs/operators/_beam_operator.py
Comment thread tests/obs/operators/test_beam_operator.py Outdated
@ArtemBasyrov
Copy link
Copy Markdown
Author

Implemented changes, suggested by the copilot review

Copy link
Copy Markdown

Copilot AI left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pull request overview

Copilot reviewed 6 out of 6 changed files in this pull request and generated 9 comments.


💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

Comment thread src/furax/math/sht.py
Comment thread src/furax/math/sht.py Outdated
Comment thread src/furax/math/sht.py
Comment on lines +152 to +188
class SHTRule(AbstractBinaryRule):
"""Algebraic rule reducing ``Map2Alm @ Alm2Map`` to an identity.

The composition *analysis after synthesis* (``Map2Alm @ Alm2Map``) is
exact when the band-limit ``lmax`` is consistent with the HEALPix
resolution ``nside``, so the rule collapses it to an
:class:`~furax.IdentityOperator` on the alm input space.

The reverse composition ``Alm2Map @ Map2Alm`` is a band-limited
projection, **not** an identity, and is therefore not reduced.
"""

left_operator_class = Map2Alm
right_operator_class = Alm2Map

def apply(
self, op1: AbstractLinearOperator, op2: AbstractLinearOperator
) -> list[AbstractLinearOperator]:
"""Reduce ``op1 @ op2`` to an identity when the pair is Map2Alm/Alm2Map.

Args:
op1: Left operator in the composition (must be :class:`Map2Alm`).
op2: Right operator in the composition (must be :class:`Alm2Map`).

Returns:
A single-element list containing an
:class:`~furax.IdentityOperator` on the alm input space.

Raises:
NoReduction: If the operator pair does not match the expected types
or if the transforms use different ``lmax`` values.
"""
if isinstance(op1, Map2Alm) and isinstance(op2, Alm2Map):
if op1.lmax != op2.lmax:
raise NoReduction
return [IdentityOperator(in_structure=op2.in_structure)]
raise NoReduction
Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The motivation behind the rule is the mathematical expectation, hence I would rather it stays

Comment thread src/furax/math/sht.py Outdated
Comment on lines +101 to +113
first_leaf = jax.tree_util.tree_leaves(x)[0]
first_leaf = jnp.atleast_2d(first_leaf) # (nfreq, npix)
nside = jhp.npix2nside(first_leaf.shape[-1])

map2alm = Map2Alm(lmax=self.lmax, nside=nside, in_structure=self.in_structure)
alm = map2alm.mv(x)

nfreq = first_leaf.shape[0]
fl = jnp.broadcast_to(jnp.atleast_2d(self.beam_fl), (nfreq, self.lmax + 1))
alm_beam = jax.tree.map(lambda alm_leaf: _apply_fl_to_alm_leaf(alm_leaf, fl), alm)

return Alm2Map(lmax=self.lmax, nside=nside, in_structure=map2alm.out_structure).mv(alm_beam)

Comment on lines +162 to +181
def mv(self, x: PyTree[Inexact[Array, 'nfreq npix']]) -> PyTree[Inexact[Array, 'nfreq npix']]:
"""Apply per-Stokes beams to input sky maps.

Args:
x: Input PyTree whose leaves have shape ``(nfreq, npix)``. The
PyTree structure must match that of ``self.beam_fl``.

Returns:
Beam-smoothed PyTree with the same structure and leaf shapes as *x*.
"""
first_leaf = jax.tree_util.tree_leaves(x)[0]
first_leaf = jnp.atleast_2d(first_leaf) # (nfreq, npix)
nside = jhp.npix2nside(first_leaf.shape[-1])

map2alm = Map2Alm(lmax=self.lmax, nside=nside, in_structure=self.in_structure)
alm = map2alm.mv(x)

alm_beam = jax.tree.map(_apply_fl_to_alm_leaf, alm, self.beam_fl)
return Alm2Map(lmax=self.lmax, nside=nside, in_structure=map2alm.out_structure).mv(alm_beam)

Comment thread src/furax/math/sht.py
Comment thread src/furax/math/sht.py Outdated
Comment thread src/furax/math/sht.py Outdated
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.

2 participants