Merge healpix mapmaker #964

earosenberg · 2024-09-24T10:59:49Z

Add healpix support to the sotodlib mapmaking module.
Includes some re-factoring of demod_mapmaker to minimize code duplication

chervias · 2024-09-26T11:43:17Z

We changed the tod_grouping file to obs_grouping. Can you fix that please?

earosenberg · 2024-09-27T13:47:19Z

We changed the tod_grouping file to obs_grouping. Can you fix that please?

I had just forgotten to delete tod_grouping I think; fixed now, thanks!

mhasself

A few comments, but then I realized this is based on the SAT mapmaking branch, not master.

Do we want this all to go into master now, or do you want to do healpix -> SAT_mapmaking first, or healpix -> master first?

sotodlib/coords/pmat.py

earosenberg · 2024-10-01T08:32:05Z

A few comments, but then I realized this is based on the SAT mapmaking branch, not master.

Do we want this all to go into master now, or do you want to do healpix -> SAT_mapmaking first, or healpix -> master first?

My thought was that this is similar in spirit to #891, so could go directly to master (and perhaps be merged from there into SAT-mapmaking). But I suppose if that is still the main branch for mapmaking the most important thing is to put it there, so we could go the other way. Happy to go with whatever you and @chervias think is easiest/best.

mhasself · 2024-10-07T17:36:51Z

@earosenberg I think it will be fine to PR this towards master rather than a mapmaking branch.

Can you rebase + squash the commit history here, please? Currently it's deeply intertwined with other branches, but then includes reversion of some changes from those branches; so it would be better to present a simpler view of what's proposed here.

chervias · 2024-10-07T17:40:15Z

I deleted the PR for the SAT_mapmaking branch, we are doing a refactor and we'll go through another route, so please merge into master.

mhasself

Thanks for pushing this forward! Lots of comments. We might need to brainstorm the right way to enrich demod_mapmaker.py classes to support healpix without too much copy paste.

We don't have enough unit testing for Pmat... that's beyond your scope, I guess -- but if you want to take a stab at inventing that, it would be welcome.

sotodlib/coords/healpix_utils.py

sotodlib/coords/pmat.py

sotodlib/mapmaking/demod_mapmaker.py

earosenberg · 2024-10-14T21:08:34Z

OK addressed above comments, did the re-factor, and added some documentation

earosenberg · 2024-10-15T16:00:37Z

Added some (simple) tests for demod_mapmaker.py and pmat.py

mhasself

Ok! I think this is looking good. A few comments about docs mostly.

Thanks for extending the existing tests -- that paves the way for more / better in the future.

Can you remind me what the use of qp_kwargs has been? It is a bit suspicious to have that, and det_weights, added to the DemodMapmaker.add_obs. You'd think that noise_model could somehow capture the effects of det_weights, and geom (or the class constructor, more generally) could handle whatever qp_kwargs is doing. I have stopped looking into this in the spirit of getting this review in in finite time, but perhaps you and @chervias could consider that, from the perspective of eventually unforking the mapmaker & signal classes.

mhasself · 2024-10-15T16:12:14Z

sotodlib/coords/healpix_utils.py

+Full map: A normal healpix map (in NEST ordering). Numpy arrays with shape (...,npix); can have any number of leading dimensions.
+Tiled map: A len(nTile) list of tiles. None for inactive tiles and a numpy array of shape (..., npix/nside) for active tiles.
+Compressed map: Numpy array of shape (nActiveTile, ..., npix/nside). Requires additional active tile information to recover a full map.
+The position of the tile axis in a compressed map can be controlled with the tile_axis_pos argument.


Thanks for the docs -- if you could "fill' these at 78 characters or so I would greatly appreciate it. It seems like you're doing one sentence per line, but rst/sphinx does not recognize that... so neither do I! Just fill them into a paragraph and let the punctuation guide the reader.

For these definitions, easiest is to just make it an rst list, with indentation to make multi-line entries:

- Full map: blah blah yada yada, then more text on next line about full map. - Tiled map: and son on. - ...

mhasself · 2024-10-15T16:14:44Z

sotodlib/coords/healpix_utils.py

+Functions are provided to transform between "full", "compressed", and "tiled" Healpix maps.
+This follows a scheme in which the full map is broken up into individual tiles defined by a lower resolution map.
+The setup ("geometry") is defined by the base nside, nside_tile of the tiling, and the ordering.
+Tiled maps must use "NEST" ordering in this scheme. If no tiling is used, "RING" is allowed.


Maybe mention that nTile = 12 * nside_tile**2... I know it's buried in the function later...

mhasself · 2024-10-15T16:21:37Z

sotodlib/coords/healpix_utils.py

+      compressed: compressed map
+      active_tile_list: len(nTile) boolean array of active tiles
+      tile_axis_pos: Index for the tile axis in the input compressed array


I think these args were copied from another function and not updated...

mhasself · 2024-10-15T16:29:11Z

sotodlib/coords/pmat.py

+      object, with attributes .shape and .wcs;
+      a pixell.tilemap.TileGeometry (if tiled);
+      or a Healpix geometry with attributes nside, nside_tile, and
+      ordering.


reformat so paragraph width is the same as the rest of this module.

mhasself · 2024-10-15T16:30:00Z

sotodlib/coords/pmat.py

-      a pixell.tilemap.TileGeometry.
+      object, with attributes .shape and .wcs;
+      a pixell.tilemap.TileGeometry (if tiled);
+      or a Healpix geometry with attributes nside, nside_tile, and


Maybe just say "or a Healpix geometry as defined in coords.healpix_utils".

mhasself · 2024-10-15T16:34:23Z

sotodlib/coords/pmat.py

+            return map
+
+    def _get_hp_tile_list(self, tiled_arr=None):
+        """Get len(nTile) boolean array of whether tiles are active. None if un-tiled.


"For healpix maps, ..."

mhasself · 2024-10-15T17:09:11Z

sotodlib/mapmaking/demod_mapmaker.py

-            self.rhs = enmap.zeros((Nsplits, ncomp)     +shape, wcs, dtype=dtype)
-            self.div = enmap.zeros((Nsplits,ncomp,ncomp)+shape, wcs, dtype=dtype)
-            self.hits= enmap.zeros((Nsplits,)+shape, wcs, dtype=dtype)
+        return cls(shape, wcs, comm, comps, name, ofmt, output, ext, dtype, sys, recenter, tile_shape, tiled, Nsplits, singlestream)


Pass parameters by name, please!

mhasself · 2024-10-15T17:09:27Z

sotodlib/mapmaking/demod_mapmaker.py

+            filter+bin mapmaking, i.e. map from obs.signal rather than from obs.dsT,
+            obs.demodQ, obs.demodU
+        """
+        return cls(None, None, None, comps, name, ofmt, output, ext, dtype, None, None, None, False, Nsplits, singlestream, nside, nside_tile)


Pass params by name.

mhasself · 2024-10-15T17:15:36Z

sotodlib/mapmaking/demod_mapmaker.py

+                    m = np.reshape(m, (np.product(m.shape[:-1]), m.shape[-1]), order='C') # Flatten wrapping axes; healpy.write_map can't handle >2d array
+                hp.write_map(oname, m.view(self.dtype), nest=(self.hp_geom.ordering=='NEST'), partial=write_partial_fits, overwrite=True)
+            elif self.ext == "npy":
+                np.save(oname, m.view(self.dtype))


I think we'd be wise to write both the array and some supporting info the projection. NEST is essential for understanding the data. NSIDE is helpful. Eventually we might throw in the celestial frame (galactic / equatorial).

The ordering and nside information is automatically included in the header for fits.
I don't think there is any way of including it for .npy files (?), except for forcibly adding it to the filename. (Aside: we can remove the npy option if it's not useful)

Like this, basically:

np.save('test.npy', {'nside': 10, 'data': np.zeros(1000)}, allow_pickle=True)

Got it, thanks for the help!

mhasself · 2024-10-15T17:15:50Z

sotodlib/mapmaking/demod_mapmaker.py

+                if h5py is None:
+                    raise ValueError("Cannot save healpix map as hdf5; h5py could not be imported. Install h5py or save as npy or fits")
+                with h5py.File(oname, 'w') as f:
+                    dset = f.create_dataset("data", m.shape, dtype=self.dtype, data=m)


Same here -- add attributes to this dataset with important projection info.

Done, added ordering and nside as attributes

earosenberg · 2024-10-16T10:09:45Z

Can you remind me what the use of qp_kwargs has been?

I am happy to leave or remove this as you and @chervias think. This was put in place to enable changing the settings of qpoint -- primarily so they could match TOAST settings. The point was to allow to_map match an input map that was scanned with TOAST at default settings. I think this is potentially useful, but maybe a sufficiently niche application that it's not worth complicating the code with.
Also, if we remove it from demod_mapmaker should I take it out of pmat as well?

You'd think that noise_model could somehow capture the effects of det_weights

I think this is correct. The main function of det_weights now is to disable inverse variance weighting for testing, but this can easily be achieved with Nmat. I'll remove it, thanks for the suggestion.

earosenberg · 2024-10-16T20:43:54Z

That should be the most recent comments done -- thanks for the continued feedback. For qp_kwargs I thought it seemed sensible to remove it from demod_mapmaker but leave it in pmat, but happy to change that if there are any objections.

By the way, I also have a unit test file for healpix_utils, mostly checking that the conversions between different healpix formats (tiled, compressed, full) are working right. Wasn't sure if you want this to be included, but will add it if so.

chervias · 2024-10-25T01:25:14Z

Besides fixing the conflict, is there anything stopping us from merging? I need this to be merged first so that I can start moving towards merging #973.

earosenberg · 2024-10-25T08:31:28Z

Not from my side -- was waiting for the go-ahead before doing a final rebase and resolving the conflict.

mhasself · 2024-10-25T15:19:23Z

Not from my side -- was waiting for the go-ahead before doing a final rebase and resolving the conflict.

This is good to go other than the conflict. Please rebase, fix conflict, squash ... in some order ... and I'll approve.

earosenberg · 2024-10-27T12:56:22Z

Rebase/conflict/squash all done.
NB the new test_pmat.py is failing post-rebase on tiled rectpix maps. It seems that the new enmapify in remove_weights and invert_weights doesn't work for the tiled case. I can remove those lines from the test if it needs to pass before merging, but it seemed like an issue to address outside this PR.

mhasself · 2024-10-28T14:27:33Z

Rebase/conflict/squash all done. NB the new test_pmat.py is failing post-rebase on tiled rectpix maps. It seems that the new enmapify in remove_weights and invert_weights doesn't work for the tiled case. I can remove those lines from the test if it needs to pass before merging, but it seemed like an issue to address outside this PR.

Ok ... ya, just suppress that test (comment it out) and we can fix it separately.

earosenberg · 2024-10-28T15:38:36Z

Ok ... ya, just suppress that test (comment it out) and we can fix it separately.

Great, done

mhasself

Quick merge it before something else breaks

earosenberg requested review from mhasself and chervias September 24, 2024 11:00

earosenberg marked this pull request as ready for review September 25, 2024 14:00

mhasself requested changes Sep 27, 2024

View reviewed changes

sotodlib/coords/pmat.py Outdated Show resolved Hide resolved

sotodlib/coords/pmat.py Outdated Show resolved Hide resolved

sotodlib/coords/pmat.py Outdated Show resolved Hide resolved

sotodlib/coords/pmat.py Outdated Show resolved Hide resolved

mhasself self-requested a review October 7, 2024 17:35

earosenberg force-pushed the merge-healpix-mapmaker branch from f20f1bd to c4ace59 Compare October 8, 2024 12:48

mhasself requested changes Oct 10, 2024

View reviewed changes

mhasself requested changes Oct 15, 2024

View reviewed changes

earosenberg force-pushed the merge-healpix-mapmaker branch from b1ce615 to b482e34 Compare October 27, 2024 12:49

mhasself mentioned this pull request Oct 28, 2024

Tiled rectpix pmat broken by recent wcs fixes #1020

Open

Add Healpix projection to coords/pmat.py and mapmaking/demod_mapmaker.py

fad90e0

earosenberg force-pushed the merge-healpix-mapmaker branch from b482e34 to fad90e0 Compare October 28, 2024 15:22

mhasself approved these changes Oct 28, 2024

View reviewed changes

earosenberg merged commit 0177b36 into master Oct 28, 2024
4 checks passed

earosenberg deleted the merge-healpix-mapmaker branch October 31, 2024 08:41

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Merge healpix mapmaker #964

Merge healpix mapmaker #964

earosenberg commented Sep 24, 2024

chervias commented Sep 26, 2024

earosenberg commented Sep 27, 2024

mhasself left a comment

earosenberg commented Oct 1, 2024

mhasself commented Oct 7, 2024

chervias commented Oct 7, 2024

mhasself left a comment

earosenberg commented Oct 14, 2024

earosenberg commented Oct 15, 2024

mhasself left a comment

mhasself Oct 15, 2024

mhasself Oct 15, 2024

mhasself Oct 15, 2024

mhasself Oct 15, 2024

mhasself Oct 15, 2024

mhasself Oct 15, 2024

mhasself Oct 15, 2024

mhasself Oct 15, 2024

mhasself Oct 15, 2024

earosenberg Oct 16, 2024 •

edited

Loading

mhasself Oct 16, 2024

earosenberg Oct 16, 2024

mhasself Oct 15, 2024

earosenberg Oct 16, 2024

earosenberg commented Oct 16, 2024 •

edited

Loading

earosenberg commented Oct 16, 2024

chervias commented Oct 25, 2024

earosenberg commented Oct 25, 2024

mhasself commented Oct 25, 2024

earosenberg commented Oct 27, 2024

mhasself commented Oct 28, 2024

earosenberg commented Oct 28, 2024

mhasself left a comment

Merge healpix mapmaker #964

Merge healpix mapmaker #964

Conversation

earosenberg commented Sep 24, 2024

chervias commented Sep 26, 2024

earosenberg commented Sep 27, 2024

mhasself left a comment

Choose a reason for hiding this comment

earosenberg commented Oct 1, 2024

mhasself commented Oct 7, 2024

chervias commented Oct 7, 2024

mhasself left a comment

Choose a reason for hiding this comment

earosenberg commented Oct 14, 2024

earosenberg commented Oct 15, 2024

mhasself left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

earosenberg Oct 16, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

earosenberg commented Oct 16, 2024 • edited Loading

earosenberg commented Oct 16, 2024

chervias commented Oct 25, 2024

earosenberg commented Oct 25, 2024

mhasself commented Oct 25, 2024

earosenberg commented Oct 27, 2024

mhasself commented Oct 28, 2024

earosenberg commented Oct 28, 2024

mhasself left a comment

Choose a reason for hiding this comment

earosenberg Oct 16, 2024 •

edited

Loading

earosenberg commented Oct 16, 2024 •

edited

Loading