more proper handling of map.pixel vs pixwin[sv] in get_pseudo2datavec.py and get_spectra_from_maps.py

zatkins2 · zatkins2 · commit 20af447465f5 · 2025-12-15T11:23:26.000-05:00
diff --git a/project/SO/pISO/python/get_mcm_bbl_and_pseudosignal.py b/project/SO/pISO/python/get_mcm_bbl_and_pseudosignal.py
@@ -46,7 +46,6 @@
     raise ValueError("Unkown 'type' value! Must be either 'Dl' or 'Cl'")
 binning_file = d["binning_file"]
 binned_mcm = d["binned_mcm"]
-mk_pseudospectra = False
 
 if d["use_toeplitz_mcm"] == True:
     assert args.old, 'can only toeplitz with pspy for now' # FIXME
diff --git a/project/SO/pISO/python/get_pseudo2datavec.py b/project/SO/pISO/python/get_pseudo2datavec.py
@@ -49,7 +49,7 @@
         # being the same survey
         templates[sv] = so_map.read_map(d[f"window_kspace_{sv}_{maps[sv][0]}"])
             
-        if d[f"pixwin_{sv}"]["pix"] == "CAR":
+        if templates[sv].pixel == "CAR":
             filter_dicts[sv] = d[f"k_filter_{sv}"]
         else:
             raise NotImplementedError('can only kspace filter CAR maps')
@@ -97,13 +97,15 @@
 
     # need to splice mbl_inv into spectra-ordered arrays
     # TODO: consider disk-space, memory (could be sparse)
+    # FIXME: function assumes TT, TE, TB, ET, BT, EE, EB, BE, BB order
     pseudo2datavec = so_mcm.get_spec2spec_array_from_spin2spin_array(mbl_inv, dense=True)
 
     # get the inv_kspace matrix for this array cross, if necessary
     if apply_kspace_filter:
         inv_kspace_mat = np.linalg.inv(kspace_transfer_matrix[spec_name]) 
 
         # apply the inv_kspace matrix to mbl_inv to get data operator
+        # FIXME: assumes inv_kspace_mat in TT, TE, TB, ET, BT, EE, EB, BE, BB order
         pseudo2datavec = inv_kspace_mat @ pseudo2datavec
 
     # get the pixwin for healpix, if necessary
diff --git a/project/SO/pISO/python/get_spectra_from_maps.py b/project/SO/pISO/python/get_spectra_from_maps.py
@@ -212,20 +212,23 @@
     # FIXME: this will not work for SO LF which has a different template despite
     # being the same survey
     templates[sv] = so_map.read_map(d[f"window_kspace_{sv}_{maps[sv][0]}"])
-        
-    if d[f"pixwin_{sv}"]["pix"] == "CAR" and deconvolve_pixwin:
-        wy, wx = enmap.calc_window(templates[sv].data.shape,
-                                   order=d[f"pixwin_{sv}"]["order"])
-        wy = wy.astype(np.float32)
-        wx = wx.astype(np.float32)
-        pixwins[sv] = (wy[:, None] * wx[None, :])
-        inv_pixwins[sv] = pixwins[sv] ** (-1)
-
-    if apply_kspace_filter:
-        filter_dicts[sv] = d[f"k_filter_{sv}"]
-        filters[sv] = kspace.get_kspace_filter(templates[sv],
-                                               filter_dicts[sv],
-                                               dtype=np.float32)
+    
+    # NOTE: a map may a CAR map but have a HEALPIX pixwin, in which case we may
+    # kspace filter it, but want to use a HEALPIX pixwin
+    if templates[sv].pixel == "CAR":
+        if d[f"pixwin_{sv}"]["pix"] == "CAR" and deconvolve_pixwin:
+            wy, wx = enmap.calc_window(templates[sv].data.shape,
+                                    order=d[f"pixwin_{sv}"]["order"])
+            wy = wy.astype(np.float32)
+            wx = wx.astype(np.float32)
+            pixwins[sv] = (wy[:, None] * wx[None, :])
+            inv_pixwins[sv] = pixwins[sv] ** (-1)
+
+        if apply_kspace_filter:
+            filter_dicts[sv] = d[f"k_filter_{sv}"]
+            filters[sv] = kspace.get_kspace_filter(templates[sv],
+                                                filter_dicts[sv],
+                                                dtype=np.float32)
 
 # get spectrum-level auxiliary data products
 spec_name_list = pspipe_list.get_spec_name_list(d, delimiter="_")
@@ -329,14 +332,14 @@
 
         window_tuple = (win_T, win_pol)
 
-        if win_T.pixel == "CAR": # FIXME: might not match d[f"pixwin_{sv}"]["pix"], there should only be one
-            if apply_kspace_filter or deconvolve_pixwin:
+        if win_T.pixel == "CAR":
+            if apply_kspace_filter or (deconvolve_pixwin and d[f"pixwin_{sv}"]["pix"] == "CAR"):
                 win_kspace = so_map.read_map(d[f"window_kspace_{sv}_{m}"])
             if apply_kspace_filter:
                 filter = filters[sv]
                 weighted_filter = filter_dicts[sv]["weighted"]
             if deconvolve_pixwin and d[f"pixwin_{sv}"]["pix"] == "CAR":
-                inv_pwin = inv_pixwins[sv] # if d[f"pixwin_{sv}"]["pix"] == "CAR" else None # FIXME: see above
+                inv_pwin = inv_pixwins[sv]
 
         cal, pol_eff = d[f"cal_{sv}_{m}"], d[f"pol_eff_{sv}_{m}"]
 
@@ -347,11 +350,7 @@
             else:
                 split_idx = k - 1 # we put signal first because it is always correlated, so has biggest memory footprint
             
-            if win_T.pixel == "CAR": # FIXME: might not match d[f"pixwin_{sv}"]["pix"], there should only be one
-
-                ###################################
-                # INJECT MAPS
-                ###################################
+            if win_T.pixel == "CAR":
 
                 # data injection
                 if which == 'data':
@@ -385,7 +384,7 @@
                         else:
                             split.write_map(f"{sim_map_dir}" + f"noise_sim_map{tag}_{sv}_{m}_set{split_idx}_{iii:05d}.fits")
 
-                if apply_kspace_filter and deconvolve_pixwin and d[f"pixwin_{sv}"]["pix"] == "CAR":
+                if apply_kspace_filter and (deconvolve_pixwin and d[f"pixwin_{sv}"]["pix"] == "CAR"):
                     if k == 0:
                         log.info(f"[Rank {so_mpi.rank}, Mapset {iii}] Apply kspace filter and inv pixwin on {sv}, {m}")
                     split = kspace.filter_map(split,
@@ -404,7 +403,7 @@
                                               weighted_filter=weighted_filter,
                                               use_ducc_rfft=True)
 
-                elif deconvolve_pixwin:
+                elif deconvolve_pixwin and d[f"pixwin_{sv}"]["pix"] == "CAR":
                     if k == 0:
                         log.info(f"[Rank {so_mpi.rank}, Mapset {iii}] WARNING: inv pixwin but no kspace filter on {sv}, {m}")
                     split = so_map.fourier_convolution(split,
@@ -416,9 +415,6 @@
                         log.info(f"[Rank {so_mpi.rank}, Mapset {iii}] WARNING: no kspace filter and no inv pixwin on {sv}, {m}")
 
             elif win_T.pixel == "HEALPIX":
-                ###################################
-                # INJECT MAPS
-                ###################################
 
                 # data injection
                 if which == 'data':
@@ -496,11 +492,14 @@
 
                 # compute specific cls with higher precision, save memory overall by
                 # doing this per spectrum. start_at_zero=False to match pspy convention
+                # TODO: test if speed penalty of alm np.complex128 conversion 
+                # is worth the memory saved (takes ~14s per spectrum)
                 _, pseudo_dict = so_spectra.get_spectra_pixell(master_alms[sv1, m1, snk1].astype(np.complex128),
                                                                master_alms[sv2, m2, snk2].astype(np.complex128),
                                                                spectra=spectra,
                                                                apply_pspy_cut=True)
                 
+                # FIXME: assumes pseudo2datavec in spectra order
                 pseudovec = so_spectra.spec_dict2vec(pseudo_dict, spectra)
                 datavec = pseudo2datavec @ pseudovec
                 
diff --git a/project/SO/pISO/python/release/add_20251201_mock_sims_to_hippo.py b/project/SO/pISO/python/release/add_20251201_mock_sims_to_hippo.py
@@ -0,0 +1,134 @@
+from os.path import join as opj
+from henry import Henry
+import hippometa as hm 
+
+client = Henry()
+sim_dir = '/scratch/gpfs/SIMONSOBS/users/zatkins/projects/lat-iso/piso_old/tiger_deep56_sim/sim_factory/sims/maps'
+
+# collection layers:
+# Mock LAT ISO-SV1 Sims
+# - Mock LAT ISO-SV1 (1000x) Sims
+#   - Mock LAT ISO-SV1 (1000x) tube freq Sims
+#     - Mock LAT Deep56 Signal Map (1000x) tube freq {(Tags)}
+#     - Mock LAT Deep56 Noise Map (1000x) tube freq (Set y)
+
+top_collection = client.new_collection(
+    name='Mock LAT ISO-SV1 Sims',
+    description='Mock simulations of the LAT-ISO SV1 data. Each simulation ' \
+    'includes the four MF tubes i1, i3, i4, and i6, and the two UHF tubes c1 ' \
+    'and i5. Each tube and frequency band includes variations of signal maps ' \
+    '(with/without random systematics, with/without map-level lensing) and ' \
+    'four-way split noise maps.\n\n' \
+    'See [this page](https://simonsobs.atlassian.net/wiki/spaces/PRO/pages/1452802049/LAT+ISO+Mock+Simulations) ' \
+    'for more technical information.',
+    products=[],
+)
+
+for sim_num in range(10):
+    sim_collection = client.new_collection(
+        name=f'Mock LAT ISO-SV1 (1000{sim_num}) Sims',
+        description='Mock simulations of the LAT-ISO SV1 data for sim index ' \
+        f'1000{sim_num}. Each simulation includes the four MF tubes i1, i3, ' \
+        'i4, and i6, and the two UHF tubes c1 and i5. Each tube and frequency ' \
+        'band includes variations of signal maps (with/without random ' \
+        'systematics, with/without map-level lensing) and four-way split noise' \
+        'maps.\n\n' \
+        'See [this page](https://simonsobs.atlassian.net/wiki/spaces/PRO/pages/1452802049/LAT+ISO+Mock+Simulations) ' \
+        'for more technical information.',
+        products=[],
+    )
+
+    for mapname in ['i1_f090', 'i1_f150', 'i3_f090', 'i3_f150', 'i4_f090', 
+                    'i4_f150', 'i6_f090', 'i6_f150', 'c1_f220', 'c1_f280', 
+                    'i5_f220', 'i5_f280']:
+        tube, freq = mapname.split('_')
+        tube_freq_collection = client.new_collection(
+            name=f'Mock LAT ISO-SV1 (1000{sim_num}) {tube} {freq} Sims',
+            description='Mock simulations of the LAT-ISO SV1 data for sim index ' \
+            f'1000{sim_num}, {tube} {freq} only. Includes variations of signal ' \
+            'maps (with/without random systematics, with/without map-level ' \
+            'lensing) and four-way split noise maps.\n\n' \
+            'See [this page](https://simonsobs.atlassian.net/wiki/spaces/PRO/pages/1452802049/LAT+ISO+Mock+Simulations) ' \
+            'for more technical information.',
+            products=[],
+        )
+
+        # signal maps
+        for tag in ['_', '_syst_', '_lens_', '_syst_lens_']:
+            map_file = f'signal_sim_map{tag}LAT_{tube}_{freq}_1000{sim_num}.fits'
+
+            if map_file == 'signal_sim_map_syst_LAT_i1_f150_10009.fits':
+                signal_map_set = client.pull_product('69387f74f22a066989fc54e3', realize_sources=False)
+            else:
+                nametag = {'_': '',
+                        '_syst_': ' (Systematics)',
+                        '_lens_': ' (Map-level Lensing)',
+                        '_syst_lens_': ' (Systematics, Map-level Lensing)'
+                        }[tag]
+                
+                desctag = {'_': '',
+                        '_syst_': ' Includes a random systematics realization.',
+                        '_lens_': ' Includes map-level lensing with lenspyx.',
+                        '_syst_lens_': ' Includes a random systematics realization and map-level lensing with lenspyx.'
+                        }[tag]
+
+                name = f'Mock LAT Deep56 Signal Map (1000{sim_num}) {tube} {freq}{nametag}'
+                desc = f'Mock signal simulation of the LAT-ISO SV1 data, for ' \
+                f'{tube} {freq}. The index of this simulation is 1000{sim_num}. ' \
+                f'Same map for each of the four-way splits.{desctag} See ' \
+                '[this page](https://simonsobs.atlassian.net/wiki/spaces/PRO/pages/1452802049/LAT+ISO+Mock+Simulations) ' \
+                'for more technical information.'
+
+                signal_map_set = client.new_product(
+                    name=name,
+                    description=desc,
+                    metadata=hm.MapSet(telescope='lat', instrument=tube,
+                                    frequency=freq, release='lat-iso-sv1',
+                                    patch='deep56',
+                                    polarization_convention='COSMO',
+                                    pixelisation='cartesian'),
+                    sources={'map': {'path': opj(sim_dir, map_file),
+                                    'description': 'the simulated map'}}
+                                    )
+
+            tube_freq_collection.append(signal_map_set)
+
+        # noise maps
+        for k in range(4):
+            map_file = f'noise_sim_map_LAT_{tube}_{freq}_set{k}_1000{sim_num}.fits'
+
+            if map_file == 'noise_sim_map_LAT_i1_f150_set0_10009.fits':
+                noise_map_set = client.pull_product('69387ef3f22a066989fc54dd', realize_sources=False)
+            else:
+                desctag = {0: ' first',
+                        1: ' second',
+                        2: ' third',
+                        3: ' fourth'
+                        }[k]
+
+                name = f'Mock LAT Deep56 Noise Map (1000{sim_num}) {tube} {freq} (Set {k})'
+                desc = f'Mock noise simulation of the LAT-ISO SV1 data, for ' \
+                f'{tube} {freq},{desctag} split of the four-way split. The index ' \
+                f'of this simulation is 1000{sim_num}. See ' \
+                '[this page](https://simonsobs.atlassian.net/wiki/spaces/PRO/pages/1452802049/LAT+ISO+Mock+Simulations) ' \
+                'for more technical information.'
+
+                noise_map_set = client.new_product(
+                    name=name,
+                    description=desc,
+                    metadata=hm.MapSet(telescope='lat', instrument=tube,
+                                    frequency=freq, release='lat-iso-sv1',
+                                    patch='deep56',
+                                    polarization_convention='COSMO',
+                                    pixelisation='cartesian', split=k),
+                    sources={'map': {'path': opj(sim_dir, map_file),
+                                    'description': 'the simulated map'}}
+                                    )
+
+            tube_freq_collection.append(noise_map_set)
+        
+        sim_collection.append(tube_freq_collection)
+    
+    top_collection.append(sim_collection)
+
+client.push(top_collection)
