Add subscan operations to preprocess #1028
Conversation
Add basic statistics, PSDs, and glitches on subscans.
earosenberg
force-pushed
the
subscan-preproc
branch
from
475c87d to
2e9b247
Compare
There was a problem hiding this comment.
Thanks @earosenberg, it'll be really nice to have this infrastructure for subscan operations. I put inline comments reflecting the conversation we had on Friday plus a few more things as I went through more closely.
sotodlib/preprocess/pcore.py
Outdated
| if (ii == 0) and isinstance(out[k], RangesMatrix): # Treat like dets | ||
| if a in full._axes: | ||
| _, fs, ns = full[a].intersection(new[a], return_slices=True) | ||
| else: | ||
| fs = range(new[a].count) | ||
| ns = range(new[a].count) | ||
| oidx.append(fs) | ||
| nidx.append(ns) | ||
| else: # Treat like samps | ||
| oidx.append(slice(None)) | ||
| nidx.append(slice(None)) |
There was a problem hiding this comment.
Can you expand the inline comment to explain what this is for/why you added it so that when something breaks and I come back here to figure out I why I remind myself why this looks so weird?
sotodlib/tod_ops/flags.py
Outdated
| subscan_aman.wrap('start_time', tt[start_inds], [(0, 'subscans')]) | ||
| stop_time = aman.timestamps[ss_ind[:-1,1]] | ||
| last_time = [2*tt[-1] - tt[-2] if end_inds[-1] == tt.size else tt[end_inds[-1]]] | ||
| subscan_aman.wrap('stop_time', np.concatenate([stop_time, last_time]), [(0, 'subscans')]) |
There was a problem hiding this comment.
It doesn't seem like you use these (start/stop_time) anywhere and I don't think we want to wrap these because in general these can change if you restrict samples which will be reflected in subscan_flags and the corresponding start/stop times can always be computed from those flags and timestamps since the both have axes aligned with samps
sotodlib/preprocess/processes.py
Outdated
| @@ -791,13 +858,20 @@ def calc_and_save(self, aman, proc_aman): | |||
| calc_aman = core.AxisManager(aman.dets, aman.samps) | |||
| calc_aman.wrap('turnarounds', ta, [(0, 'dets'), (1, 'samps')]) | |||
|
|
|||
| self.save(proc_aman, calc_aman) | |||
| if ('merge_subscans' not in self.calc_cfgs) or (self.calc_cfgs['merge_subscans']): | |||
| subscan_aman = aman.subscan_info | |||
There was a problem hiding this comment.
can this just go in calc_aman? It seems like all other methods look for subscan_info inside of aman not proc_aman since FlagTurnarounds.process has a call to get_turnaround_flags that repopulates aman.subscan_info.
sotodlib/preprocess/processes.py
Outdated
| name = "stats" | ||
| def __init__(self, step_cfgs): | ||
| self.signal = step_cfgs.get('signal', 'signal') | ||
| self.wrap = step_cfgs.get('wrap', 'stats') |
There was a problem hiding this comment.
can we make the default to be tod_stats or something similar since we already have hwpss_stats and azss_stats
sotodlib/preprocess/processes.py
Outdated
| calc: | ||
| stat_names: ["median", "std"] | ||
| split_subscans: False # optional | ||
| mask: # optional, for cutting a power spectrum in frequency |
There was a problem hiding this comment.
how about we call this psd_mask and rename TODStats just to GetStats or something similar that doesn't imply this only operates on the TOD.
sotodlib/tod_ops/fft_ops.py
Outdated
| @@ -353,6 +390,8 @@ def fit_noise_model( | |||
| f_max=100, | |||
| merge_name="noise_fit_stats", | |||
| merge_psd=True, | |||
| freq_spacing=None, | |||
| approx_fit=False, | |||
There was a problem hiding this comment.
Let's get rid of approx_fit and just have the behavior of fit in the preprocess method do what approx_fit is doing there. This is essentially doing what's in calc_wn without the fknee and alpha estimates so perhaps we expand calc_wn to include those, call that fn to set p0 and when you want "approx_fit" you just call that fn.
There was a problem hiding this comment.
Sounds good.
It's a little confusing that calc_wn returns the square root of median(psd) whereas fit_noise_model doesn't take the root; do we want to change that? Also seems cleaner to add a new function estimate_fknee but I can mash it all into calc_wn if you prefer
There was a problem hiding this comment.
Yes, this bothers me quite a lot that we (probably I) decided to change the convention between calc_wn and fit_noise_model. I think that they should both be in squared units (like fit_noise_model is now). I see what you're saying about wanting to have another function like estimate_fknee but I think my preference is to have one function called estimate_noise_model which returns exactly the same shaped output as fit_noise_model i.e. an axis shaped (dets x noise_coeffs). If we want to keep calc_wn fn alive for backwards compatibility we could just have it call estimate_noise_model (so that we don't have to update definitions of median/mean/etc in multiple places) and just return the white noise coefficients but I would keep everything in the squared units for consistency.
Another related issue was brought up by @ykyohei in #995 where we really should be always taking the PSDs w/ overlap = 0, storing the nperseg in the calc_psd output so that we can have the calc_wn/any fn that takes a median over psd bins correct for this "median bias".
It's maybe ok to address that all here but also ok opening up a new PR for these fixes and reviewing that separately, I'll leave it up to you.
There was a problem hiding this comment.
After some consideration I think estimate_noise_model either needs to be more robust then just the initial guess of the fit_noise_model or at least have some error quantification. Otherwise you can have catastrophic mis-fits -- this is also a problem with fit_noise_model but I assume it could be caught with the covariance. I've just deleted approx_fit for now, will open a new PR (or at least an issue) to deal with these other problems.
sotodlib/tod_ops/flags.py
Outdated
| """ | ||
| Returns an axis manager with information about subscans. | ||
| This includes direction, start time, stop time, and a ranges matrix (subscans samps) | ||
| True inside each subscan. Subscans are defined excluding turnarounds. |
There was a problem hiding this comment.
Let's follow the docstring format adding in
Arguments
-----------
arg: type
description
Returns
--------
return_variable: type
description
sotodlib/tod_ops/flags.py
Outdated
| return arr[ndslice] | ||
|
|
||
| def wrap_info(aman, info_aman_name, info, info_names, merge=True): | ||
| """Wrap info into an aman. Assumed to be (dets,) or (dets, subscans)""" |
There was a problem hiding this comment.
This function at first glance seems to be quite general but I think it's somewhat specific to wrapping in the stats info appropriately choosing whether or not to include subscan axes. Can you just call it wrap_stats or wrap_and_check_subscans or something similarly specific and give the doc string a bit more description.
sotodlib/tod_ops/flags.py
Outdated
| stat_names : list | ||
| List of strings identifying which statistics to run. | ||
| split_subscans : bool | ||
| (Optional). If True statistics will be computed on subscans. Assumes aman.subscans exists already. |
| signal = aman[signal] | ||
| if split_subscans: | ||
| if mask is not None: | ||
| raise ValueError("Cannot mask samples and split subscans") |
There was a problem hiding this comment.
It should be possible to do a union of the subscan flags + a secondary flag if they're both ranges() aligned along samps. Ok to leave for later but I'd just create an issue for it to track it.
There was a problem hiding this comment.
Added support for discontinuous ranges to apply_rng but didn't implement the full masking. Added #1046 .
|
Addressed all above points except where noted. Also added new features:
|
Add basic statistics, PSDs, and glitches on subscans.