Speed up the binner, compute the HWP angle once, fix nside_centering and singular pixels in the destriper

[paganol]

Two small speedups found while profiling a 3-year, 75 Hz run on Galileo100, where scan_map and the binner each took ~1600 s, mostly computing pointings on the fly.

Binner: analytic inversion of the 3×3 matrices. _solve_binning called np.linalg.cond, solve and inv for every pixel. Now np.linalg.eigvalsh only computes the eigenvalues (eigvalsh), and the inverse uses the adjugate.

nside 512, 3.1M pixels	before	after
_solve_binning	1.9 s	0.35 s

On 200k test pixels (also scaled by 1e-12 and 1e8), the pixels marked as UNSEEN are identical, and the solutions differ by cond·ε: 8e-14 for cond < 1e3, 3e-10 up to 1e6. I first tried closed formulas for the eigenvalues (estimate_cond_number, used by the destriper), but they lose too much precision on the smallest eigenvalue of singular matrices (1–2 hits per pixel), and accepted pixels that np.linalg.eigvalsh rejects.

Pointings: the HWP angle once per observation. PointingProvider.get_pointings recomputed the HWP angle at every call, although it does not depend on the detector; the map-makers, scan_map, etc. call it once per detector. The new parameter compute_hwp_angle (default True) lets _get_pointings_array, Observation.get_pointings("all") and get_hitmap skip it.

8 detectors × 1 day, 75 Hz	before	after
per-detector pointing loop	4.19 s	3.16 s

New test test_smart_pointings_skip_hwp_angle.

Bug fix: nside_centering with a single detector. Observation.get_pointings(det_idx, nside_centering=...) with an integer index returned the raw pointings: the centering was applied only for a list or "all". The map-makers and the beam convolution were not affected, as they center the pointings in _get_pointings_array. The centering is now in a small helper used by both branches; the new test test_smart_pointings_nside_centering fails without the fix.

Bug fix: singular pixels in the destriper. The destriper used estimate_cond_number to decide which pixels are valid. With weights of order 1, it accepts ~25% of singular matrices (1–2 hits). With realistic weights (1/(f_s·NET²), NET in µK√s) it rejects all of them, but only thanks to an absolute threshold (1e-10) on the smallest eigenvalue, so the result depended on the units of the NETs. The binner and the destriper now share _is_well_conditioned, which uses eigvalsh and is scale-independent. The destriper reference tests are unchanged. The check costs 0.5 s instead of 0.03 s at nside 512, once per destriper run. New test test_is_well_conditioned.

Abort the MPI job on an unhandled exception. If one process raises, the others wait for it in their next MPI call (e.g., the collective free of the shared-memory windows at exit), and the job hangs until the wall-time limit (a known mpi4py issue). As suggested by @anand-avinash, litebird_sim/mpi.py now installs the exception hook used in BrahMap when MPI is enabled with more than one process: it prints the traceback and calls MPI_Abort, like python -m mpi4py. mpi.rst documents it. I checked locally that a 2-rank job where one rank raises now terminates in ~1 s instead of hanging. Test: test_mpi_excepthook.

The full suite passes locally.

🤖 Generated with Claude Code

[github-actions]

Coverage report

Click to see where and how coverage changed

File	Statements	Missing	Coverage	Coverage (new stmts)	Lines missing
litebird_sim
mpi.py					61-63, 76-79, 83
observation_utilities.py
observations.py					1142
pointings.py
litebird_sim/mapmaking
binner.py					88-119
common.py					379, 393-394
destriper.py					422-429
Project Total

This report was generated by python-coverage-comment-action

[mj-gomes]

Hi @paganol, just a small comment which does not interfere in any way with the functionality of the changes.

You mention LAPACK here and on the code comments, but since LAPACK is a fortran package and people reading LBS code will not necessarily be aware of it, I think we should stick to calling things in their "pythonic way" (numpy's linalg). Of course linalg is based on LAPACK and BLAS, but, at least for me, it was confusing since I did not see any reference to LAPACK elsewhere in LBS.

[anand-avinash]

...so if one process fails, the job hangs until the wall-time limit

This is a known issue of mpi4py and may happen irrespective of shared-memory usage. mpi4py indeed recommends using -m mpi4py option while running python scripts.

Another solution is using an exception hook. It allows clean exception handling while running mpi jobs without having to use -m mpi4py option. It calls COMM_WORLD.Abort(1) from the problematic process to kill the mpi job elegantly. I have used it in BrahMap, which you can check for the reference.