Add support for geom objects to be rendered with galactic or ecliptic frames by erykoff · Pull Request #185 · LSSTDESC/healsparse

erykoff · 2026-10-23T22:59:52Z

This PR also adds more convenient ways to clear pixels, including del m[geom] and some optimizations made possible when using pixel ranges.

Also, thanks to some tips from @sidneymau on some memory inefficiencies he has seen I did some investigation into the memory usage when combining very large packed boolean maps and found:

There was a problem in the way I was reading in packed boolean arrays from fits which meant that the arrays weren't resizable without a temporary (slow) memory doubling. This was fixed.
When doing invert there were unnecessary temporaries being created that could be fixed by doing inversion in place.
The bit counting was refactored to take advantage of numpy 2 features which makes it much faster and use much less memory. So doing n_valid or get_valid_area() on a bit-packed map is 12x faster for DP2-sized maps.
numpy 2 is now required.

… not append.

kadrlica · 2026-10-24T15:19:20Z

I don't know the motivation for del m[geom]. My personal take is:

del m[geom] seems overboard to me.
m[geom] = None and m &= geom seem reasonable
m.clear_pixels(geom) is nice as an explicit interface

This allows faster memory resizing when combining maps.

…ory. With claude assistance to save memory on collapsing the last axis.

sidneymau

A few minor comments, but everything looks reasonable to me overall!

sidneymau · 2026-10-28T04:25:11Z

    m += healsparse.Circle(ra=202.0, dec=21.0, radius=0.5, value=10)


+Using Geometric Shapes To Set And Erase Pixels


From the documentation, it's not clear what is meant by "erasing" a pixel -- does that mean to set the value to UNSEEN? To free the associated memory?

I will rewrite the documentation to make this clear. And I should say "clear" pixels, not erase, because that's what the code is called. And what it does is (a) set all pixels to the sentinel value (which is typically UNSEEN for floating point maps, and False for boolean maps, and either 0 or -MAXINT for integer maps) (b) never expand the coverage map (because all of those are already invalid/sentinel values). The associated memory is not freed unless you defragment the map. I'll add all this to the docs.

sidneymau · 2026-10-28T04:29:02Z

+        in this frame, and output will be equatorial/celestial coordinates.
    """
-    def __init__(self, *, ra, dec, radius, value, nside_render=None):
+    def __init__(self, *, ra, dec, radius, value, nside_render=None, geom_frame="equatorial"):


I would think frame would make more sense than geom_frame throughout -- it's implicit that the frame is the local frame for this object already

Makes sense! I'll update.

sidneymau · 2026-10-28T04:32:44Z

        self._n_valid = None

+        # Convert to an array, no-copy if possible.
+        _pixels = np.atleast_1d(np.asarray(pixels))


I think atleast_1d already converts objects to appropriate array-like objects, so asarray is probably unnecessary (and redundant) here. Did you find this was necessary?

I think you are right.

sidneymau · 2026-10-28T04:36:45Z

+        start, stop = pixel_ranges[:, 0], pixel_ranges[:, 1]
+        cov0 = np.right_shift(start, shift)
+        cov1 = np.right_shift(stop - 1, shift)  # inclusive; no overrun at the top
+        nseg = cov1 - cov0 + 1


Is this guaranteed to be nonnegative? I can't seem to break that, so it probably is

As long as the pixel ranges are in healpix range ... but now I'm not sure if that's explicitly checked or if it'll just raise an error when called into hpgeom.

sidneymau · 2026-10-28T04:38:01Z

+        cov0 = np.right_shift(start, shift)
+        cov1 = np.right_shift(stop - 1, shift)  # inclusive; no overrun at the top
+        nseg = cov1 - cov0 + 1
+        rng = np.repeat(np.arange(len(start)), nseg)


Using rng as the variable name for anything other than a random number generator is bound to be confusing in the future

sidneymau · 2026-10-28T04:41:23Z

                             (key.__class__))

+    def clear_pixels(self, key):
+        """Clear part of a healsparse map.


Similar to my comment about "Erase" earlier -- what is actually happening here? m[key] = None is also shorthand for something more specific. Is the memory associated with each of these keys freed? Or are the values densely set to UNSEEN?

I will expand the docstring along with the documentation.

sidneymau · 2026-10-28T04:49:19Z

                summand += np.sum(last_unpacked[0], dtype=np.int64)
            if mid_data is not None:
-                summand += np.sum(self._bit_count(mid_data), dtype=np.int64)
+                summand += np.sum(self._last_axis_bit_counts(mid_data, mid_data.size * 8))


Should this use a dtype itemsize instead of 8?

Ah, so the packed boolean array must always beuint8. That gives the most flexibility and efficiency for packing bits. I will add a self._itemsize or something just so it is always clear.

But the changes here on the PR involve the "count number of bits set", not the "check the bits that are set" (which I haven't noticed a problem with but I'll profile that a bit). When I originally wrote this code numpy2 was only in beta, and therefore I couldn't use the new numpy2 bitwise_count routine. This is much faster than the LUT method I was doing before. And it can be made even FASTER if you use wide words (simply viewing the uint8 array as 64-bit words) which then has two effects: (a) the temporary output is 8x smaller (less memory overhead and smaller allocations); (b) numpy will use a special optimized bit-counting instruction on processors that support it. I actually think that the biggest win is moving to the numpy bitwise_count (which is in c rather than via python); second is the memory reduction; third is the special instruction. But all that is somewhat academic, why not do all the optimizations available in numpy 2?

sidneymau · 2026-10-28T04:52:56Z

@@ -681,16 +683,31 @@ def _test_bits_at_locs(self, locs):

        return self._data[_locs // 8] & (1 << (_locs % 8).astype(np.uint8)) != 0


another place where np.dtype(np.int64).itemsize might be useful for generality? E.g., if we ever want a super high resolution map that needs int128 pixels, or if there is ever a need for optimizing smaller maps with smaller bit-packed integers? Or maybe that is overkill...

erykoff added 11 commits October 22, 2026 17:01

Fix crash if a single bare pixel was sent to update_values_pix().

e9c4f19

Add test for setting a single pixel to update_values_pix().

9d1e863

Add geom support for __setitem__().

9297b2f

Add tests for geom in __setitem__.

29ddec9

Add __delitem__ and clear_pixels for easier pixel clears.

4dee946

Add tests for __delitem__ and clear_pixels().

f4d499e

Update update_values_pix so that a boolean map and-ed with False does…

9ac5f7f

… not append.

Optimize generating a map from a geom object.

c509cff

Add frame rotation support to geometry library.

332ad06

Add tests for frame rotations for geometry objects.

10d52b1

Add documentation on geometric frames and deleting geometric objects.

769860e

erykoff added 7 commits October 25, 2026 12:09

Remove superfluous __delitem__ for clearing pixels.

00ad830

Refactor _update_values_pixel_ranges() with help of claude.

a8990f4

Fix initialization of GeomUnion and GeomIntersection.

2ecec41

Add tests for other GeomUnion and GeomIntersection initialization.

5b6defc

Add functionality to crop maps to geom objects efficiently.

e6ceffc

Add tests for cropping functionality.

9ad2cab

Add documentation for crop method.

c53a059

erykoff force-pushed the geom_frames branch 2 times, most recently from bbeafcb to 43823f4 Compare October 27, 2026 18:07

erykoff added 6 commits October 27, 2026 13:10

Use resize instead of ravel when reading bit-packed maps.

2c57790

This allows faster memory resizing when combining maps.

Remove duplicate call to check coverage pixels.

84bea94

Update minimum numpy to 2.0

f22350f

Refactor bit counting to take advantage of numpy 2.0 and conserve mem…

63eec3c

…ory. With claude assistance to save memory on collapsing the last axis.

Make boolean inversions much more memory efficient.

5ba2aef

Use ravel for healpix read.

60d6ba8

erykoff force-pushed the geom_frames branch from baee14b to 60d6ba8 Compare October 27, 2026 20:41

sidneymau approved these changes Oct 28, 2026

View reviewed changes

Rename geom_frame to simply frame.

30892c4

erykoff added 2 commits October 28, 2026 09:22

Clean up some variable names.

f685cbe

Update docstrings/documentation on what it means to "clear" pixels.

81b788a

erykoff merged commit 258972e into main Oct 28, 2026
5 checks passed

erykoff deleted the geom_frames branch October 28, 2026 19:43

		m += healsparse.Circle(ra=202.0, dec=21.0, radius=0.5, value=10)


		Using Geometric Shapes To Set And Erase Pixels

		@@ -681,16 +683,31 @@ def _test_bits_at_locs(self, locs):

		return self._data[_locs // 8] & (1 << (_locs % 8).astype(np.uint8)) != 0

Conversation

erykoff commented Oct 23, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

kadrlica commented Oct 24, 2026

Uh oh!

sidneymau left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

erykoff commented Oct 23, 2026 •

edited

Loading