Skip to content

Modify mainly source_flags so that we can use multiple sources and taua. #1235

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 1 commit into from
Jun 5, 2025

Conversation

yoshinori-0778
Copy link
Contributor

This PR is for pipeline planet mapmaking, which contains four changes:

  • source_flags including planet.py:
    • These change allows us to use multiple sources and taua like the block below
    • I also modify source selection so that we can cut detectors that is not observing the source. And I add options for cutting detectors based on e.g. the Moon.
    • I add function for making Inverce mask of source mask.
     - name : "source_flags"
       calc:
         mask: {'shape': 'circle',
               'xyr': [0, 0, 2.0]}
         center_on: ['jupiter', 'tauA'] # list of str
         res: 2 # arcmin
         max_pix: 4000000 # max number of allowed pixels in map
         distance: 0 # max distance of footprint from source in degrees
         inv_flag: True
       save: True
       select: # optional
         kind: "any"
         select_source: ['jupiter']
    
  • trim_flag_edge: I add this new function for triming TOD because I sometimes have a problem when I use filter_for_source (pca filter). This problem is due to interplay between gapfill and flags. If we have flags of True at the end of TOD, gapfill try to fill those based on only the prior TOD whereas gapfill generally try to fill based on both prior and posterior TODs. This can degrade the quality of the gapfill. This function trim both the beginning and end of the TOD so that gapfill is never based on only one side.
  • glitch: This change allows us to set any flags name.
  • correct_iir_params: fix a bug

I checked that these works as expected, but I am happy to modify them more if needed.

@yoshinori-0778 yoshinori-0778 requested a review from msilvafe May 28, 2025 19:37
@yoshinori-0778
Copy link
Contributor Author

Kyohei pointed me out that

  • glitch modification is not needed because original code already has. I did not noticed that, so I will take it back after confirmation.
  • Also it might be better to use Jxxx+yyy for taua. I will consider better way to use it.
    It is not needed to review until I finish the updates above.

@@ -119,7 +119,11 @@ def get_scan_q(tod, planet, boresight_offset=None, refq=None):
el = np.median(tod.boresight.el[::10])
az = np.median(tod.boresight.az[::10])
t = (tod.timestamps[0] + tod.timestamps[-1]) / 2
if isinstance(planet, str):
if isinstance(planet, (list, tuple)):
planet_name, ra, dec = planet
Copy link
Contributor

Choose a reason for hiding this comment

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

planet_name unused replace with _

Copy link
Contributor Author

Choose a reason for hiding this comment

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

fixed this.

Copy link
Contributor

Choose a reason for hiding this comment

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

Not a big deal but it doesn't actually look like you changed planet_name -> _

@@ -171,7 +171,6 @@ class GlitchDetection(_FracFlaggedMixIn, _Preprocess):
Example configuration block::

- name: "glitches"
glitch_name: "my_glitches"
Copy link
Contributor

Choose a reason for hiding this comment

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

Ya this has been updated, I think you must have been testing off an old version of the code.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I did test with an old version. I took back this modify.

Comment on lines 1258 to 1259
if self.calc_cfgs.get('inv_flag') is not None:
if self.calc_cfgs['inv_flag']:
Copy link
Contributor

Choose a reason for hiding this comment

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

Should just be able to do: if self.calc_cfgs.get('inv_flag') and not 2 layers of if-statement. If the inv_flag arg isn't present it will return None and None is treated the same as False.

Similarly don't need is not None in L1283 and 1294

Copy link
Contributor Author

Choose a reason for hiding this comment

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

That is good to know. I modified these.

Comment on lines 1220 to 1230
else:
updated_list = []
for isource in SOURCE_LIST:
if isinstance(isource, str):
if isource in source_list:
updated_list.append(isource)
elif isinstance(isource, tuple):
isource_name = isource[0]
if isource_name in source_list:
updated_list.append(isource)
source_list = updated_list
Copy link
Contributor

Choose a reason for hiding this comment

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

If sourcelist and SOURCE_LIST were the same type (i.e. all str or all tuple) then you could do np.intersect1d(source_list, SOURCE_LIST) instead of these ~10 lines right? I don't know all of the history behind why there's two possible types but perhaps you can make a helper function in coords that convers the list into the same format as SOURCE_LIST so that you can call:

source_list = coord.planets.normalize_source_types(source_list)
source_list = np.intersect1d(source_list, SOURCE_LIST)

Or just wrap it into one function inside like: coords.planets.get_source_list_overlap(source_list) and internally it can handle the types and intersection with SOURCE_LIST. Probably could write a similar helper for L1237-1245. Ideally we keep the stuff in preprocess/processes slim (and I could imagine this is a thing you may wanna do outside of a preprocess run). Although I'm not very well versed in using the coords.planets module especially for these extra added sources specified as tuples so if this doesn't make sense you can ignore.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I made new function get_source_list_fromstr in planets.py and change scripts so that new function is used.


self.save(proc_aman, source_aman)

def save(self, proc_aman, source_aman):
if self.save_cfgs is None:
return
if self.save_cfgs:
proc_aman.wrap("source_flags", source_aman)
if self.calc_cfgs.get('source_flags_name') is not None:
proc_aman.wrap(self.calc_cfgs.get('source_flags_name', None), source_aman)
Copy link
Contributor

Choose a reason for hiding this comment

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

Like other functions let's introduce source_flags_name at the top level of the class and add an _init_ function then call self.source_flags_name here. I know we use calc_cfgs in select and process in a few other of the preprocess methods but we would like to avoid doing that where unneccesary.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Modified these based on your comments. This is also good to know, thank you!


source_list = np.atleast_1d(self.calc_cfgs.get('center_on', 'planet'))
if source_list == ['planet']:
select_list = np.atleast_1d(self.select_cfgs["select_source"])
Copy link
Contributor

Choose a reason for hiding this comment

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

Give an example of this in the example config entry(ies) in the docstring. Currently I think if you passed in self.select_cfgs = True as is done in the current docstring example this will fail bc you're missing select_source for that reason to observe backwards compatibility choose a reasonable default--probably all source flags.

elif self.select_cfgs["kind"] == "cut":
keep = ~has_any_cuts(source_flags[source])
else:
raise ValueError(f"Unknown kind of selection: {self.select_cfgs['kind']}")
meta.restrict("dets", meta.dets.vals[keep])
Copy link
Contributor

Choose a reason for hiding this comment

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

Do one restrict at the end. Start with an empty keep flag then keep += each of these over the select_list then run the restrict on that total keep at the end. Also I think a useful extension here would be not just to cut based on any or all criteria but if > X% of samples are flagged for a detector.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Modified these based on your comments. I made new function has_ratio_cuts and flag_cut_select. default values of new parameters written in yaml file are set for backward compatibility.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

@mmccrackan I noticed that recent PR #1116 modified these part and conflict with my version now. how should I do? This is the first time to have a conflict...
Should I revert my change and make another PR for this?

Copy link
Contributor

Choose a reason for hiding this comment

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

I think it should be fine to add a small change here. Basically you'll just want to pass in_place=True into select(). If True, return the restricted meta. Otherwise, return the final keep after accumulating all the flags you want in it.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Thank you for your comment, I slightly modified my modification to incorporate your modification. Also rebased to avoid conflict.

Comment on lines 2068 to 2234
class TrimFlagEdge(_Preprocess):
"""Trim edge until given flags of all detectors are False
To find first and last sample id that has False (i.e., no flags applied) for all detectors.
This is for avoiding glitchfill problem for data whose edge has flags of True.
"""
name = 'trim_flag_edge'

def process(self, aman, proc_aman, sim=False):
flags = aman.flags.get(self.process_cfgs.get('flags'))
trimst = np.where(~np.any(flags.mask(), axis = 0))[0][0]
trimen = np.where(~np.any(flags.mask(), axis = 0))[0][-1]
aman.restrict('samps', (aman.samps.offset + trimst,
aman.samps.offset + trimen))
proc_aman.restrict('samps', (proc_aman.samps.offset + trimst,
proc_aman.samps.offset + trimen))
Copy link
Contributor

Choose a reason for hiding this comment

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

We handle this in the glitch detection step by not allowing samples to be flagged within some buffer number of samples from the beginning or end of the TOD could that be an extension to this instead? Also it should be able to compute this from the ranges tuples instead of completely unpacking from the RangesMatrices into boolean masks twice (that conversion can be slow always faster if you can work within Ranges/sparse arrays). Also need to add an example config entry for this.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Modified these so that we calculate these indices with range. New version is much faster than old version by about order of three.
Also very good to know how you handle this in the glitch detection. I agree that we can do the same treatment in the filter_for_sources but I would like to add this class because it need additional study. I could try this type of treatment later.

@msilvafe
Copy link
Contributor

Also @davidvng and @mmccrackan have contributed to this part of the repo before perhaps they want to take a look before we merge.

@mmccrackan mmccrackan self-requested a review May 30, 2025 14:21
@yoshinori-0778
Copy link
Contributor Author

Thank you for a lot of comments. I have fixed codes as you suggested except for the one for triming. I am considering how I should modify it. Once I finish this, I will push all changes.

@yoshinori-0778
Copy link
Contributor Author

I modified scritps based on Max's comments. Thank you Max.

@yoshinori-0778 yoshinori-0778 force-pushed the planet_mapmaker_common branch from 78ef202 to 1838e78 Compare June 2, 2025 04:33
@yoshinori-0778
Copy link
Contributor Author

I rebased and slightly modified it to avoid conflict and incorporate @mmccrackan's modification #1116.

@yoshinori-0778
Copy link
Contributor Author

@msilvafe @mmccrackan This is ready for review again. Would you mind if you do review when you have a time?

@msilvafe
Copy link
Contributor

msilvafe commented Jun 4, 2025

@msilvafe @mmccrackan This is ready for review again. Would you mind if you do review when you have a time?

Sounds good, I'll take a look this morning. For the future, if you click the "re-request review" button on the top right of the PR I'll get notified that you're ready for me to look over it again.

@yoshinori-0778
Copy link
Contributor Author

I see, I did not know that button, good to know, thank you.
Thank you for taking time for review as well!

Copy link
Contributor

@msilvafe msilvafe left a comment

Choose a reason for hiding this comment

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

A couple inline comments but I'm approving the changes. Thanks for making all of the requested updates. I think we will benefit from having the updates you added to flagman and more examples of manipulating RangesMatrices without expanding to boolean arrays.

@@ -119,7 +119,11 @@ def get_scan_q(tod, planet, boresight_offset=None, refq=None):
el = np.median(tod.boresight.el[::10])
az = np.median(tod.boresight.az[::10])
t = (tod.timestamps[0] + tod.timestamps[-1]) / 2
if isinstance(planet, str):
if isinstance(planet, (list, tuple)):
planet_name, ra, dec = planet
Copy link
Contributor

Choose a reason for hiding this comment

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

Not a big deal but it doesn't actually look like you changed planet_name -> _

considered to have too many cuts.
"""

if not isinstance(ratio, (float, int)):
Copy link
Contributor

Choose a reason for hiding this comment

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

how is it interpreted if it's an int? I.e. if I pass in 10 should it be 10% or is this basically to capture 0 or 1?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I am assuming this is basically to capture 0 or 1. I modified this if part so that an error show up when ratio < 0 or ratio > 1 is passed.

Returns:
min_idx, max_idx: minmum and maximum indices that has False flag across all detectros.
"""
max_val = max(arr.complement().ranges()[:, 1].max() for arr in flags)
Copy link
Contributor

Choose a reason for hiding this comment

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

No need to change anything here I just found this pretty cool. Had never seen this syntax before where you can pass the generator directly into the max function. Happy to learn some new python from this review :)


Args:
flags (RangesMatrix): An instance of so3g.proj.RangesMatrix indicating flagged time ranges.
Returns:
Copy link
Contributor

Choose a reason for hiding this comment

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

Since flag_cut_select can flip the logic of the flag do you want this function to also be able to look for edges with either logic?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I do not have an situation that we want to find edges that all have True now, so I would like to leave this for now. But thank you for the suggestion.

len_samps = flag.shape[-1]
return np.array([np.sum(np.diff(x.ranges()))/len_samps > ratio for x in flag])

def flag_cut_select(flags, cut, kind):
Copy link
Contributor

Choose a reason for hiding this comment

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

I think a better name for cut would be invert since this is used in source flags and the standard flag definition across sotodlib is samples to be excluded = True, samples to be kept = False so invert would flip that logic. Also I think cut (or now invert) should be an optional argument that defaults to the standard logic we use across the library.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Sure, I changed the name to invert and modified corresponding parts.

@yoshinori-0778 yoshinori-0778 force-pushed the planet_mapmaker_common branch from 1838e78 to 0180e30 Compare June 4, 2025 20:20
@yoshinori-0778
Copy link
Contributor Author

Thank you for the review.
I pushed new version that is reflected your comment and squash all commits because there was another conflict.
I am going to merge this end of today.

@yoshinori-0778 yoshinori-0778 merged commit 21202c1 into master Jun 5, 2025
5 checks passed
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.

3 participants