From 0175b8e96533700dfab96ffce93cd4f920492d8f Mon Sep 17 00:00:00 2001 From: adrien-laposta Date: Fri, 10 Oct 2025 05:33:24 -0400 Subject: [PATCH 01/12] compute and use data from different pipeline stages in sim filtering --- sotodlib/preprocess/pcore.py | 30 ++++++- sotodlib/preprocess/preprocess_util.py | 59 ++++++++----- sotodlib/preprocess/processes.py | 116 +++++++++++++++++++------ 3 files changed, 152 insertions(+), 53 deletions(-) diff --git a/sotodlib/preprocess/pcore.py b/sotodlib/preprocess/pcore.py index 8443e170f..12bff2292 100644 --- a/sotodlib/preprocess/pcore.py +++ b/sotodlib/preprocess/pcore.py @@ -38,7 +38,8 @@ def __init__(self, step_cfgs): self.select_cfgs = step_cfgs.get("select") self.plot_cfgs = step_cfgs.get("plot") self.skip_on_sim = step_cfgs.get("skip_on_sim", False) - def process(self, aman, proc_aman, sim=False): + self.use_data_aman = step_cfgs.get("use_data_aman", False) + def process(self, aman, proc_aman, sim=False, data_aman=None): """ This function makes changes to the time ordered data AxisManager. Ex: calibrating or detrending the timestreams. This function will use any configuration information under the ``process`` key of the @@ -56,6 +57,9 @@ def process(self, aman, proc_aman, sim=False): sim: Bool False by default when analyzing data. Should be True when doing Transfer Function simulations and determining which steps should be run. + data_aman: AxisManager (Optional) + An AxisManager containing the preprocessed data to be used by + this process. """ if self.process_cfgs is None: return aman, proc_aman @@ -435,8 +439,14 @@ def extend(self, index, other): super().extend( [self._check_item(item) for item in other]) def __setitem__(self, index, item): super().__setitem__(index, self._check_item(item)) + def __getitem__(self, index): + result = super().__getitem__(index) + if isinstance(index, slice): + return Pipeline(result) + else: + return result - def run(self, aman, proc_aman=None, select=True, sim=False, update_plot=False): + def run(self, aman, proc_aman=None, select=True, sim=False, update_plot=False, data_amans=None): """ The main workhorse function for the pipeline class. This function takes an AxisManager TOD and successively runs the pipeline of preprocessing @@ -472,6 +482,11 @@ def run(self, aman, proc_aman=None, select=True, sim=False, update_plot=False): given ``proc_aman`` is ``aman.preprocess``. This assumes ``process.calc_and_save()`` has been run on this aman before and has injested flags and other information into ``proc_aman``. + data_amans: dict (Optional) + A dictionary of AxisManagers with keys (step, process.name) + filled with AxisManager processed up to step-1. This is used + to pre-load all data AxisManager which could be required when + processing simulations (e.g. to provide a T2P template) Returns ------- @@ -523,7 +538,16 @@ def run(self, aman, proc_aman=None, select=True, sim=False, update_plot=False): if sim and process.skip_on_sim: continue self.logger.debug(f"Running {process.name}") - aman, proc_aman = process.process(aman, proc_aman, sim) + if (data_amans is not None) and process.use_data_aman: + try: + data_aman = data_amans[step, process.name] + except KeyError: + raise KeyError(f"Requested to use data AxisManager for process {process.name} but not found in data_amans") + else: + if process.use_data_aman and sim: + raise ValueError(f"Process {process.name} requested to use data_aman but none was provided to Pipeline.run()") + data_aman = None + process.process(aman, proc_aman, sim, data_aman) if run_calc: aman, proc_aman = process.calc_and_save(aman, proc_aman) process.plot(aman, proc_aman, filename=os.path.join(self.plot_dir, '{ctime}/{obsid}', f'{step+1}_{{name}}.png')) diff --git a/sotodlib/preprocess/preprocess_util.py b/sotodlib/preprocess/preprocess_util.py index 11d83539c..cc9993d01 100644 --- a/sotodlib/preprocess/preprocess_util.py +++ b/sotodlib/preprocess/preprocess_util.py @@ -390,7 +390,8 @@ def load_and_preprocess(obs_id, configs, context=None, dets=None, meta=None, def multilayer_load_and_preprocess(obs_id, configs_init, configs_proc, dets=None, meta=None, no_signal=None, - logger=None, init_only=False): + logger=None, init_only=False, + stop_for_sims=False): """Loads the saved information from the preprocessing pipeline from a reference and a dependent database and runs the processing section of the pipeline for each. @@ -421,6 +422,11 @@ def multilayer_load_and_preprocess(obs_id, configs_init, configs_proc, Optional. Logger object or None will generate a new one. init_only: bool Optional. If True, do not run the dependent pipeline. + stop_for_sims: bool + Optinal. If True, will stop before each step of the pipeline + with the flag `use_data_aman` set to True. The intended use is + to prepare all necessary data products that cannot be stored in + the preprocessing database, to process simulations. """ if logger is None: @@ -482,9 +488,29 @@ def multilayer_load_and_preprocess(obs_id, configs_init, configs_proc, if 'valid_data' in aman.preprocess: aman.preprocess.move('valid_data', None) aman.preprocess.merge(proc_aman.preprocess) - pipe_proc.run(aman, aman.preprocess, select=False) + if stop_for_sims: + batch_idx = [ + (step, process.name) + for step, process in enumerate(pipe_proc) + if process.use_data_aman + ] + batch_idx = [(0, pipe_proc[0].name)] + batch_idx + pipes = {} + for idx in range(len(batch_idx)-1): + start, start_name = batch_idx[idx] + end, end_name = batch_idx[idx+1] + pipes[end, end_name] = pipe_proc[start:end] + + out_amans = {} + loc_aman = aman.copy() + for (step, name), pipe in pipes.items(): + pipe.run(loc_aman, aman.preprocess, select=False) + out_amans[step, name] = loc_aman.copy() + return out_amans - return aman + else: + pipe_proc.run(aman, aman.preprocess, select=False) + return aman else: raise ValueError('Dependency check between configs failed.') @@ -492,7 +518,7 @@ def multilayer_load_and_preprocess(obs_id, configs_init, configs_proc, def multilayer_load_and_preprocess_sim(obs_id, configs_init, configs_proc, sim_map, meta=None, logger=None, init_only=False, - t2ptemplate_aman=None): + data_amans=None): """Loads the saved information from the preprocessing pipeline from a reference and a dependent database, loads the signal from a (simulated) map into the AxisManager and runs the processing section of the pipeline @@ -524,9 +550,11 @@ def multilayer_load_and_preprocess_sim(obs_id, configs_init, configs_proc, Optional. Logger object or None will generate a new one. init_only: bool Optional. Whether or not to run the dependent pipeline. - t2ptemplate_aman: AxisManager - Optional. AxisManager to use as a template for t2p leakage - deprojection. + data_amans: dict (Optional) + A dictionary of AxisManagers with keys (step, process.name) + filled with AxisManager processed up to step-1. This is used + to pre-load all data AxisManager which could be required when + processing simulations (e.g. to provide a T2P template) """ if logger is None: logger = init_logger("preprocess") @@ -580,28 +608,13 @@ def multilayer_load_and_preprocess_sim(obs_id, configs_init, configs_proc, if init_only: return aman - - if t2ptemplate_aman is not None: - # Replace Q,U with simulated timestreams - t2ptemplate_aman.wrap("demodQ", aman.demodQ, [(0, 'dets'), (1, 'samps')], overwrite=True) - t2ptemplate_aman.wrap("demodU", aman.demodU, [(0, 'dets'), (1, 'samps')], overwrite=True) - - t2p_aman = t2pleakage.get_t2p_coeffs( - t2ptemplate_aman, - merge_stats=False - ) - t2pleakage.subtract_t2p( - aman, - t2p_aman, - T_signal=t2ptemplate_aman.dsT - ) logger.info("Running dependent pipeline") proc_aman = context_proc.get_meta(obs_id, meta=aman) if 'valid_data' in aman.preprocess: aman.preprocess.move('valid_data', None) aman.preprocess.merge(proc_aman.preprocess) - pipe_proc.run(aman, aman.preprocess, sim=True) + pipe_proc.run(aman, aman.preprocess, sim=True, data_amans=data_amans) return aman else: diff --git a/sotodlib/preprocess/processes.py b/sotodlib/preprocess/processes.py index a53a1a4d7..1c1707cf6 100644 --- a/sotodlib/preprocess/processes.py +++ b/sotodlib/preprocess/processes.py @@ -30,7 +30,9 @@ class FFTTrim(_Preprocess): .. autofunction:: sotodlib.tod_ops.fft_trim """ name = "fft_trim" - def process(self, aman, proc_aman, sim=False): + def process(self, aman, proc_aman, sim=False, data_aman=None): + if data_aman is not None: + raise NotImplementedError("No support for using data AxisManager in process") start_stop = tod_ops.fft_trim(aman, **self.process_cfgs) proc_aman.restrict(self.process_cfgs.get('axis', 'samps'), (start_stop)) @@ -48,7 +50,9 @@ def __init__(self, step_cfgs): super().__init__(step_cfgs) - def process(self, aman, proc_aman, sim=False): + def process(self, aman, proc_aman, sim=False, data_aman=None): + if data_aman is not None: + raise NotImplementedError("No support for using data AxisManager in process") tod_ops.detrend_tod(aman, signal_name=self.signal, **self.process_cfgs) return aman, proc_aman @@ -275,7 +279,9 @@ def __init__(self, step_cfgs): super().__init__(step_cfgs) - def process(self, aman, proc_aman, sim=False): + def process(self, aman, proc_aman, sim=False, data_aman=None): + if data_aman is not None: + raise NotImplementedError("No support for using data AxisManager in process") field = self.process_cfgs['jumps_aman'] aman[self.signal] = tod_ops.jumps.jumpfix_subtract_heights( aman[self.signal], proc_aman[field].jump_flag.mask(), @@ -749,7 +755,9 @@ def __init__(self, step_cfgs): super().__init__(step_cfgs) - def process(self, aman, proc_aman, sim=False): + def process(self, aman, proc_aman, sim=False, data_aman=None): + if data_aman is not None: + raise NotImplementedError("No support for using data AxisManager in process") if self.process_cfgs["kind"] == "single_value": if self.process_cfgs.get("divide", False): aman[self.signal] /= self.process_cfgs["val"] @@ -926,7 +934,9 @@ def __init__(self, step_cfgs): super().__init__(step_cfgs) - def process(self, aman, proc_aman, sim=False): + def process(self, aman, proc_aman, sim=False, data_aman=None): + if data_aman is not None: + raise NotImplementedError("No support for using data AxisManager in process") if not(proc_aman[self.hwpss_stats] is None): modes = [int(m[1:]) for m in proc_aman[self.hwpss_stats].modes.vals[::2]] if sim: @@ -1002,7 +1012,9 @@ class Apodize(_Preprocess): """ name = "apodize" - def process(self, aman, proc_aman, sim=False): + def process(self, aman, proc_aman, sim=False, data_aman=None): + if data_aman is not None: + raise NotImplementedError("No support for using data AxisManager in process") tod_ops.apodize.apodize_cosine(aman, **self.process_cfgs) return aman, proc_aman @@ -1036,7 +1048,9 @@ class Demodulate(_Preprocess): name = "demodulate" - def process(self, aman, proc_aman, sim=False): + def process(self, aman, proc_aman, sim=False, data_aman=None): + if data_aman is not None: + raise NotImplementedError("No support for using data AxisManager in process") hwp.demod_tod(aman, **self.process_cfgs["demod_cfgs"]) if self.process_cfgs.get("trim_samps"): trim = self.process_cfgs["trim_samps"] @@ -1134,7 +1148,9 @@ def save(self, proc_aman, azss_stats): if self.save_cfgs: proc_aman.wrap(self.calc_cfgs["azss_stats_name"], azss_stats) - def process(self, aman, proc_aman, sim=False): + def process(self, aman, proc_aman, sim=False, data_aman=None): + if data_aman is not None: + raise NotImplementedError("No support for using data AxisManager in process") if 'subtract_in_place' in self.calc_cfgs: raise ValueError('calc_cfgs.subtract_in_place is not allowed use process_cfgs.subtract') if self.process_cfgs is None: @@ -1183,7 +1199,9 @@ class SubtractAzSSTemplate(_Preprocess): """ name = "subtract_azss_template" - def process(self, aman, proc_aman, sim=False): + def process(self, aman, proc_aman, sim=False, data_aman=None): + if data_aman is not None: + raise NotImplementedError("No support for using data AxisManager in process") process_cfgs = copy.deepcopy(self.process_cfgs) if sim: process_cfgs["azss"] = proc_aman.get(process_cfgs["azss"]) @@ -1219,7 +1237,9 @@ def __init__(self, step_cfgs): super().__init__(step_cfgs) - def process(self, aman, proc_aman, sim=False): + def process(self, aman, proc_aman, sim=False, data_aman=None): + if data_aman is not None: + raise NotImplementedError("No support for using data AxisManager in process") if self.flags is not None: glitch_flags=proc_aman.get(self.flags) tod_ops.gapfill.fill_glitches( @@ -1273,7 +1293,9 @@ def save(self, proc_aman, turn_aman): if self.save_cfgs: proc_aman.wrap("turnaround_flags", turn_aman) - def process(self, aman, proc_aman, sim=False): + def process(self, aman, proc_aman, sim=False, data_aman=None): + if data_aman is not None: + raise NotImplementedError("No support for using data AxisManager in process") tod_ops.flags.get_turnaround_flags(aman, **self.process_cfgs) return aman, proc_aman @@ -1285,7 +1307,9 @@ class SubPolyf(_Preprocess): """ name = 'sub_polyf' - def process(self, aman, proc_aman, sim=False): + def process(self, aman, proc_aman, sim=False, data_aman=None): + if data_aman is not None: + raise NotImplementedError("No support for using data AxisManager in process") tod_ops.sub_polyf.subscan_polyfilter(aman, **self.process_cfgs) return aman, proc_aman @@ -1641,7 +1665,9 @@ class HWPAngleModel(_Preprocess): """ name = "hwp_angle_model" - def process(self, aman, proc_aman, sim=False): + def process(self, aman, proc_aman, sim=False, data_aman=None): + if data_aman is not None: + raise NotImplementedError("No support for using data AxisManager in process") if (not 'hwp_angle' in aman._fields) and ('hwp_angle' in proc_aman._fields): aman.wrap('hwp_angle', proc_aman['hwp_angle']['hwp_angle'], [(0, 'samps')]) @@ -1725,7 +1751,9 @@ def __init__(self, step_cfgs): super().__init__(step_cfgs) - def process(self, aman, proc_aman, sim=False): + def process(self, aman, proc_aman, sim=False, data_aman=None): + if data_aman is not None: + raise NotImplementedError("No support for using data AxisManager in process") field = self.process_cfgs.get("noise_fit_array", None) if field: noise_fit = proc_aman[field] @@ -1951,7 +1979,9 @@ def __init__(self, step_cfgs): super().__init__(step_cfgs) - def process(self, aman, proc_aman, sim=False): + def process(self, aman, proc_aman, sim=False, data_aman=None): + if data_aman is not None: + raise NotImplementedError("No support for using data AxisManager in process") n_modes = self.process_cfgs.get('n_modes') signal = aman.get(self.signal) if aman.dets.count < n_modes: @@ -2015,7 +2045,9 @@ def __init__(self, step_cfgs): super().__init__(step_cfgs) - def process(self, aman, proc_aman, sim=False): + def process(self, aman, proc_aman, sim=False, data_aman=None): + if data_aman is not None: + raise NotImplementedError("No support for using data AxisManager in process") n_modes = self.process_cfgs.get('n_modes') signal = aman.get(self.signal) flags = aman.flags.get(self.process_cfgs.get('source_flags')) @@ -2171,9 +2203,23 @@ class SubtractT2P(_Preprocess): """ name = "subtract_t2p" - def process(self, aman, proc_aman, sim=False): - tod_ops.t2pleakage.subtract_t2p(aman, proc_aman['t2p'], - **self.process_cfgs) + def process(self, aman, proc_aman, sim=False, data_aman=None): + if data_aman is not None: + data_aman.wrap("demodQ", aman.demodQ, [(0, 'dets'), (1, 'samps')], overwrite=True) + data_aman.wrap("demodU", aman.demodU, [(0, 'dets'), (1, 'samps')], overwrite=True) + + t2p_aman = tod_ops.t2pleakage.get_t2p_coeffs( + data_aman, + merge_stats=False + ) + tod_ops.t2pleakage.subtract_t2p( + aman, + t2p_aman, + T_signal=data_aman.dsT + ) + else: + tod_ops.t2pleakage.subtract_t2p(aman, proc_aman['t2p'], + **self.process_cfgs) return aman, proc_aman class SplitFlags(_Preprocess): @@ -2230,7 +2276,9 @@ class UnionFlags(_Preprocess): """ name = "union_flags" - def process(self, aman, proc_aman, sim=False): + def process(self, aman, proc_aman, sim=False, data_aman=None): + if data_aman is not None: + raise NotImplementedError("No support for using data AxisManager in process") warnings.warn("UnionFlags function is deprecated and only kept to allow loading of old process archives. Use generalized method CombineFlags") from so3g.proj import RangesMatrix @@ -2266,7 +2314,9 @@ class CombineFlags(_Preprocess): """ name = "combine_flags" - def process(self, aman, proc_aman, sim=False): + def process(self, aman, proc_aman, sim=False, data_aman=None): + if data_aman is not None: + raise NotImplementedError("No support for using data AxisManager in process") from so3g.proj import RangesMatrix if isinstance(self.process_cfgs['method'], list): if len(self.process_cfgs['flag_labels']) != len(self.process_cfgs['method']): @@ -2315,7 +2365,9 @@ class RotateFocalPlane(_Preprocess): """ name = "rotate_focal_plane" - def process(self, aman, proc_aman, sim=False): + def process(self, aman, proc_aman, sim=False, data_aman=None): + if data_aman is not None: + raise NotImplementedError("No support for using data AxisManager in process") from sotodlib.coords import demod demod.rotate_focal_plane(aman, **self.process_cfgs) return aman, proc_aman @@ -2335,7 +2387,9 @@ class RotateQU(_Preprocess): """ name = "rotate_qu" - def process(self, aman, proc_aman, sim=False): + def process(self, aman, proc_aman, sim=False, data_aman=None): + if data_aman is not None: + raise NotImplementedError("No support for using data AxisManager in process") from sotodlib.coords import demod demod.rotate_demodQU(aman, **self.process_cfgs) return aman, proc_aman @@ -2373,7 +2427,9 @@ def save(self, proc_aman, aman): if self.save_cfgs: proc_aman.wrap('qu_common_mode_coeffs', aman['qu_common_mode_coeffs']) - def process(self, aman, proc_aman, sim=False): + def process(self, aman, proc_aman, sim=False, data_aman=None): + if data_aman is not None: + raise NotImplementedError("No support for using data AxisManager in process") if 'qu_common_mode_coeffs' in proc_aman: tod_ops.deproject.subtract_qu_common_mode(aman, self.signal_name_Q, self.signal_name_U, coeff_aman=proc_aman['qu_common_mode_coeffs'], @@ -2440,7 +2496,9 @@ class PointingModel(_Preprocess): """ name = "pointing_model" - def process(self, aman, proc_aman, sim=False): + def process(self, aman, proc_aman, sim=False, data_aman=None): + if data_aman is not None: + raise NotImplementedError("No support for using data AxisManager in process") from sotodlib.coords import pointing_model if self.process_cfgs: pointing_model.apply_pointing_model(aman) @@ -2519,7 +2577,9 @@ class CorrectIIRParams(_Preprocess): """ name = "correct_iir_params" - def process(self, aman, proc_aman, sim=False): + def process(self, aman, proc_aman, sim=False, data_aman=None): + if data_aman is not None: + raise NotImplementedError("No support for using data AxisManager in process") from sotodlib.obs_ops import correct_iir_params correct_iir_params(aman) @@ -2552,7 +2612,9 @@ class TrimFlagEdge(_Preprocess): """ name = 'trim_flag_edge' - def process(self, aman, proc_aman, sim=False): + def process(self, aman, proc_aman, sim=False, data_aman=None): + if data_aman is not None: + raise NotImplementedError("No support for using data AxisManager in process") flags = aman.flags.get(self.process_cfgs.get('flags')) trimst, trimen = core.flagman.find_common_edge_idx(flags) aman.restrict('samps', (aman.samps.offset + trimst, From d30594137e9e988685a19d6c96537d46e53fb2d6 Mon Sep 17 00:00:00 2001 From: adrien-laposta Date: Fri, 10 Oct 2025 05:58:35 -0400 Subject: [PATCH 02/12] make `skip_on_sim` a required field --- sotodlib/preprocess/pcore.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/sotodlib/preprocess/pcore.py b/sotodlib/preprocess/pcore.py index 12bff2292..37e4a374d 100644 --- a/sotodlib/preprocess/pcore.py +++ b/sotodlib/preprocess/pcore.py @@ -37,7 +37,9 @@ def __init__(self, step_cfgs): self.save_cfgs = step_cfgs.get("save") self.select_cfgs = step_cfgs.get("select") self.plot_cfgs = step_cfgs.get("plot") - self.skip_on_sim = step_cfgs.get("skip_on_sim", False) + self.skip_on_sim = step_cfgs.get("skip_on_sim") + if self.skip_on_sim is None: + raise ValueError(f"Process {step_cfgs.get('name')} must have `skip_on_sim` key set to True or False") self.use_data_aman = step_cfgs.get("use_data_aman", False) def process(self, aman, proc_aman, sim=False, data_aman=None): """ This function makes changes to the time ordered data AxisManager. From 4b9ac4a3964fb7875ae8d9e76575552ff7708570 Mon Sep 17 00:00:00 2001 From: adrien-laposta Date: Fri, 10 Oct 2025 06:14:05 -0400 Subject: [PATCH 03/12] move skip_on_sim warning to run() --- sotodlib/preprocess/pcore.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sotodlib/preprocess/pcore.py b/sotodlib/preprocess/pcore.py index 37e4a374d..0149ffbc2 100644 --- a/sotodlib/preprocess/pcore.py +++ b/sotodlib/preprocess/pcore.py @@ -38,8 +38,6 @@ def __init__(self, step_cfgs): self.select_cfgs = step_cfgs.get("select") self.plot_cfgs = step_cfgs.get("plot") self.skip_on_sim = step_cfgs.get("skip_on_sim") - if self.skip_on_sim is None: - raise ValueError(f"Process {step_cfgs.get('name')} must have `skip_on_sim` key set to True or False") self.use_data_aman = step_cfgs.get("use_data_aman", False) def process(self, aman, proc_aman, sim=False, data_aman=None): """ This function makes changes to the time ordered data AxisManager. @@ -537,6 +535,8 @@ def run(self, aman, proc_aman=None, select=True, sim=False, update_plot=False, d success = 'end' for step, process in enumerate(self): + if sim and (process.skip_on_sim is None): + raise ValueError(f"Process {process.name} missing required field `skip_on_sim`") if sim and process.skip_on_sim: continue self.logger.debug(f"Running {process.name}") From 61e1e2a6c686dc6e5768517e68e348d43bd9f1dc Mon Sep 17 00:00:00 2001 From: adrien-laposta Date: Thu, 27 Nov 2025 04:51:25 -0500 Subject: [PATCH 04/12] warning when saving multiple instances of axisman --- sotodlib/preprocess/preprocess_util.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/sotodlib/preprocess/preprocess_util.py b/sotodlib/preprocess/preprocess_util.py index cc9993d01..e591ad3ec 100644 --- a/sotodlib/preprocess/preprocess_util.py +++ b/sotodlib/preprocess/preprocess_util.py @@ -438,6 +438,21 @@ def multilayer_load_and_preprocess(obs_id, configs_init, configs_proc, configs_proc, context_proc = get_preprocess_context(configs_proc) meta_proc = context_proc.get_meta(obs_id, dets=dets, meta=meta) + # Count number of stops + if stop_for_sims: + num_stops = 0 + for process in configs_init["process_pipe"]: + if process.get("use_data_aman", False): + num_stops += 1 + for process in configs_proc["process_pipe"]: + if process.get("use_data_aman", False): + num_stops += 1 + logger.warning( + "Currently running with `stop_for_sims=True`. " + f"It will generate {num_stops} additional copies " + "of the data AxisManager with a higher memory usage." + ) + group_by_init, groups_init, error_init = get_groups(obs_id, configs_init, context_init) group_by_proc, groups_proc, error_proc = get_groups(obs_id, configs_proc, context_proc) From 5fca3c890beef6bea03013cdaa97fde441cff9a5 Mon Sep 17 00:00:00 2001 From: adrien-laposta Date: Thu, 27 Nov 2025 08:55:35 -0500 Subject: [PATCH 05/12] extend to init pipe --- sotodlib/preprocess/preprocess_util.py | 89 ++++++++++++++++++++------ 1 file changed, 69 insertions(+), 20 deletions(-) diff --git a/sotodlib/preprocess/preprocess_util.py b/sotodlib/preprocess/preprocess_util.py index e591ad3ec..a9ab1cc37 100644 --- a/sotodlib/preprocess/preprocess_util.py +++ b/sotodlib/preprocess/preprocess_util.py @@ -493,34 +493,37 @@ def multilayer_load_and_preprocess(obs_id, configs_init, configs_proc, aman = context_init.get_obs(meta_init, no_signal=no_signal) logger.info("Running initial pipeline") - pipe_init.run(aman, aman.preprocess, select=False) - if init_only: - return aman + if stop_for_sims: + out_amans_init = run_pipeline_stepgroups( + pipe_init, + aman, + run_last_step=not(init_only) + ) + if init_only: + return out_amans_init + else: + pipe_init.run(aman, aman.preprocess, select=False) + if init_only: + return aman logger.info("Running dependent pipeline") + if stop_for_sims: + aman = out_amans_init[(len(pipe_init), 'last')] proc_aman = context_proc.get_meta(obs_id, meta=aman) if 'valid_data' in aman.preprocess: aman.preprocess.move('valid_data', None) aman.preprocess.merge(proc_aman.preprocess) if stop_for_sims: - batch_idx = [ - (step, process.name) - for step, process in enumerate(pipe_proc) - if process.use_data_aman - ] - batch_idx = [(0, pipe_proc[0].name)] + batch_idx - pipes = {} - for idx in range(len(batch_idx)-1): - start, start_name = batch_idx[idx] - end, end_name = batch_idx[idx+1] - pipes[end, end_name] = pipe_proc[start:end] - - out_amans = {} - loc_aman = aman.copy() - for (step, name), pipe in pipes.items(): - pipe.run(loc_aman, aman.preprocess, select=False) - out_amans[step, name] = loc_aman.copy() + out_amans = run_pipeline_stepgroups( + pipe_proc, + out_amans_init[(len(pipe_init), 'last')], + ) + out_amans.update({ + (step, name): out_amans_init[(step, name)] + for (step, name) in out_amans_init + if name != 'last' + }) return out_amans else: @@ -530,6 +533,52 @@ def multilayer_load_and_preprocess(obs_id, configs_init, configs_proc, raise ValueError('Dependency check between configs failed.') +def run_pipeline_stepgroups(pipe, aman, run_last_step=False): + """ + Run a Pipeline object, grouping steps based on + the flag `use_data_aman` in the configuration + file. + Arguments + ---------- + pipe : Pipeline + Pipeline object to run. + aman : AxisManager + AxisManager to process. + run_last_step : bool + If True, will create a dict item containing the + AxisManager after run the full pipeline. + """ + batch_idx = [ + (step, process.name) + for step, process in enumerate(pipe) + if process.use_data_aman + ] + if batch_idx or run_last_step: + batch_idx = [(0, pipe[0].name)] + batch_idx + if run_last_step: + batch_idx += [(len(pipe), 'last')] + pipes = {} + for idx in range(len(batch_idx)-1): + start, start_name = batch_idx[idx] + end, end_name = batch_idx[idx+1] + # If asked to stop at the first process + # one needs to save the current state of + # the AxisManager + if end == 0: + pipes[end, end_name] = None + else: + pipes[end, end_name] = pipe[start:end] + out_amans = {} + loc_aman = aman.copy() + for (step, name), pipe in pipes.items(): + if pipe is not None: + pipe.run(loc_aman, aman.preprocess, select=False) + out_amans[step, name] = loc_aman.copy() + return out_amans + else: + return {} + + def multilayer_load_and_preprocess_sim(obs_id, configs_init, configs_proc, sim_map, meta=None, logger=None, init_only=False, From e29df6f362b00cc3fbc7f85040bd359d78c3b02e Mon Sep 17 00:00:00 2001 From: Adrien La Posta Date: Thu, 19 Feb 2026 14:32:09 +0000 Subject: [PATCH 06/12] adapt to fft t2p coeff fit --- sotodlib/preprocess/processes.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/sotodlib/preprocess/processes.py b/sotodlib/preprocess/processes.py index 6c6d53bfa..f6e66e488 100644 --- a/sotodlib/preprocess/processes.py +++ b/sotodlib/preprocess/processes.py @@ -2580,14 +2580,20 @@ def __init__(self, step_cfgs): super().__init__(step_cfgs) def process(self, aman, proc_aman, sim=False, data_aman=None): - if data_aman is not None: + if sim and data_aman is not None: data_aman.wrap("demodQ", aman.demodQ, [(0, 'dets'), (1, 'samps')], overwrite=True) data_aman.wrap("demodU", aman.demodU, [(0, 'dets'), (1, 'samps')], overwrite=True) - t2p_aman = tod_ops.t2pleakage.get_t2p_coeffs( - data_aman, - merge_stats=False - ) + if self.process_cfgs.get("fit_in_freq"): + t2p_aman = tod_ops.t2pleakage.get_t2p_coeffs_in_freq( + data_aman, + merge_stats=False + ) + else: + t2p_aman = tod_ops.t2pleakage.get_t2p_coeffs( + data_aman, + merge_stats=False + ) tod_ops.t2pleakage.subtract_t2p( aman, t2p_aman, From d6183d74fe6f287b30a7ac47f40f82d4dffa3c16 Mon Sep 17 00:00:00 2001 From: Adrien La Posta Date: Thu, 19 Feb 2026 14:42:53 +0000 Subject: [PATCH 07/12] add data_aman arg to processes --- sotodlib/preprocess/processes.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/sotodlib/preprocess/processes.py b/sotodlib/preprocess/processes.py index f6e66e488..e60b7724c 100644 --- a/sotodlib/preprocess/processes.py +++ b/sotodlib/preprocess/processes.py @@ -427,7 +427,9 @@ def __init__(self, step_cfgs): super().__init__(step_cfgs) - def process(self, aman, proc_aman, sim=False): + def process(self, aman, proc_aman, sim=False, data_aman=None): + if data_aman is not None: + raise NotImplementedError("No support for using data AxisManager in process") full_output = self.process_cfgs.get('full_output') if full_output: freqs, Pxx, nseg = tod_ops.fft_ops.calc_psd(aman, signal=aman[self.signal], @@ -2867,7 +2869,9 @@ def __init__(self, step_cfgs): super().__init__(step_cfgs) - def process(self, aman, proc_aman, sim=False): + def process(self, aman, proc_aman, sim=False, data_aman=None): + if data_aman is not None: + raise NotImplementedError("No support for using data AxisManager in process") scan_freq = tod_ops.utils.get_scan_freq(aman) hpf_cfg = {'type': 'sine2', 'cutoff': scan_freq, 'trans_width': scan_freq/10} filt = tod_ops.get_hpf(hpf_cfg) @@ -3231,7 +3235,9 @@ def __init__(self, step_cfgs): super().__init__(step_cfgs) - def process(self, aman, proc_aman, sim=False): + def process(self, aman, proc_aman, sim=False, data_aman=None): + if data_aman is not None: + raise NotImplementedError("No support for using data AxisManager in process") aman.move(**self.process_cfgs) return aman, proc_aman From e23958880a15ad305afd7b3d4d8f3033a3350173 Mon Sep 17 00:00:00 2001 From: Adrien La Posta Date: Fri, 20 Feb 2026 15:49:34 +0000 Subject: [PATCH 08/12] fix inconsistent ndets --- sotodlib/preprocess/preprocess_util.py | 18 ++++++++++++++---- sotodlib/preprocess/processes.py | 4 +++- 2 files changed, 17 insertions(+), 5 deletions(-) diff --git a/sotodlib/preprocess/preprocess_util.py b/sotodlib/preprocess/preprocess_util.py index 11cdeb637..eb5fb3722 100644 --- a/sotodlib/preprocess/preprocess_util.py +++ b/sotodlib/preprocess/preprocess_util.py @@ -632,7 +632,9 @@ def multilayer_load_and_preprocess(obs_id, configs_init, configs_proc, return None else: pipe_init = Pipeline(configs_init["process_pipe"], logger=logger) - aman_cfgs_ref = get_pcfg_check_aman(pipe_init) + + if not ignore_cfg_check: + aman_cfgs_ref = get_pcfg_check_aman(pipe_init) if ( ignore_cfg_check or @@ -747,6 +749,7 @@ def run_pipeline_stepgroups(pipe, aman, run_last_step=False): def multilayer_load_and_preprocess_sim(obs_id, configs_init, configs_proc, sim_map, meta=None, logger=None, init_only=False, + ignore_cfg_check=False, data_amans=None): """Loads the saved information from the preprocessing pipeline from a reference and a dependent database, loads the signal from a (simulated) @@ -779,6 +782,9 @@ def multilayer_load_and_preprocess_sim(obs_id, configs_init, configs_proc, Optional. Logger object or None will generate a new one. init_only : bool Optional. Whether or not to run the dependent pipeline. + ignore_cfg_check : bool + If True, do not attempt to validate that configs_init is the same as + the config used to create the existing init db. data_amans: dict (Optional) A dictionary of AxisManagers with keys (step, process.name) filled with AxisManager processed up to step-1. This is used @@ -811,10 +817,14 @@ def multilayer_load_and_preprocess_sim(obs_id, configs_init, configs_proc, return None else: pipe_init = Pipeline(configs_init["process_pipe"], logger=logger) - aman_cfgs_ref = get_pcfg_check_aman(pipe_init) - if check_cfg_match(aman_cfgs_ref, meta_proc.preprocess['pcfg_ref'], - logger=logger): + if not ignore_cfg_check: + aman_cfgs_ref = get_pcfg_check_aman(pipe_init) + + if ignore_cfg_check or check_cfg_match( + aman_cfgs_ref, + meta_proc.preprocess['pcfg_ref'], + logger=logger): pipe_proc = Pipeline(configs_proc["process_pipe"], logger=logger) logger.info("Restricting detectors on all proc pipeline processes") diff --git a/sotodlib/preprocess/processes.py b/sotodlib/preprocess/processes.py index e60b7724c..d5343f788 100644 --- a/sotodlib/preprocess/processes.py +++ b/sotodlib/preprocess/processes.py @@ -57,7 +57,7 @@ def __init__(self, step_cfgs): self.save_name = None super().__init__(step_cfgs) - + def process(self, aman, proc_aman, sim=False, data_aman=None): if data_aman is not None: raise NotImplementedError("No support for using data AxisManager in process") @@ -2583,6 +2583,8 @@ def __init__(self, step_cfgs): def process(self, aman, proc_aman, sim=False, data_aman=None): if sim and data_aman is not None: + + data_aman.restrict("dets", aman.dets.vals) data_aman.wrap("demodQ", aman.demodQ, [(0, 'dets'), (1, 'samps')], overwrite=True) data_aman.wrap("demodU", aman.demodU, [(0, 'dets'), (1, 'samps')], overwrite=True) From d6ef3542e5f516e531ca75886706813ee24d5e36 Mon Sep 17 00:00:00 2001 From: Adrien La Posta Date: Tue, 24 Mar 2026 10:49:44 +0000 Subject: [PATCH 09/12] ML fit for T2P coeffs with sims --- sotodlib/preprocess/processes.py | 3 +- sotodlib/tod_ops/t2pleakage.py | 144 +++++++++++++++++-------------- 2 files changed, 83 insertions(+), 64 deletions(-) diff --git a/sotodlib/preprocess/processes.py b/sotodlib/preprocess/processes.py index d5343f788..af423ec48 100644 --- a/sotodlib/preprocess/processes.py +++ b/sotodlib/preprocess/processes.py @@ -2591,7 +2591,8 @@ def process(self, aman, proc_aman, sim=False, data_aman=None): if self.process_cfgs.get("fit_in_freq"): t2p_aman = tod_ops.t2pleakage.get_t2p_coeffs_in_freq( data_aman, - merge_stats=False + merge_stats=False, + ML_fit=True ) else: t2p_aman = tod_ops.t2pleakage.get_t2p_coeffs( diff --git a/sotodlib/tod_ops/t2pleakage.py b/sotodlib/tod_ops/t2pleakage.py index 6a3804aa8..3e31d6914 100644 --- a/sotodlib/tod_ops/t2pleakage.py +++ b/sotodlib/tod_ops/t2pleakage.py @@ -304,6 +304,7 @@ def get_t2p_coeffs(aman, T_sig_name='dsT', Q_sig_name='demodQ', U_sig_name='demo def get_t2p_coeffs_in_freq(aman, T_sig_name='dsT', Q_sig_name='demodQ', U_sig_name='demodU', fs=None, fit_freq_range=(0.01, 0.1), wn_freq_range=(0.2, 1.9), subtract_sig=False, merge_stats=True, t2p_stats_name='t2p_stats', + ML_fit=False ): """ Compute the leakage coefficients from temperature (T) to polarization (Q and U) in Fourier @@ -333,6 +334,9 @@ def get_t2p_coeffs_in_freq(aman, T_sig_name='dsT', Q_sig_name='demodQ', U_sig_na Whether to merge the calculated statistics back into `aman`. Default is True. t2p_stats_name : str Name under which to wrap the output AxisManager containing statistics. Default is 't2p_stats'. + ML_fit : bool + Whether to use the ML solution (analytic) to fit for coefficients. Default is False, + which uses ODR to perform a linear fit. Returns ------- @@ -347,72 +351,86 @@ def get_t2p_coeffs_in_freq(aman, T_sig_name='dsT', Q_sig_name='demodQ', U_sig_na I_fs = rfft(aman[T_sig_name], axis=1) Q_fs = rfft(aman[Q_sig_name], axis=1) U_fs = rfft(aman[U_sig_name], axis=1) - - coeffsQ = np.zeros(aman.dets.count) - errorsQ = np.zeros(aman.dets.count) - redchi2sQ = np.zeros(aman.dets.count) - coeffsU = np.zeros(aman.dets.count) - errorsU = np.zeros(aman.dets.count) - redchi2sU = np.zeros(aman.dets.count) fit_mask = (fit_freq_range[0] < freqs) & (freqs < fit_freq_range[1]) - wn_mask = (wn_freq_range[0] < freqs) & (freqs < wn_freq_range[1]) - def leakage_model(B, x): - return B[0] * x - - model = Model(leakage_model) - - for i in range(aman.dets.count): - # fit Q - Q_wnl = np.nanmean(np.abs(Q_fs[i][wn_mask])) - x = np.real(I_fs[i])[fit_mask] - y = np.real(Q_fs[i])[fit_mask] - sx = Q_wnl / np.sqrt(2) * np.ones_like(x) - sy = Q_wnl * np.ones_like(y) - try: - data = RealData(x=x, - y=y, - sx=sx, - sy=sy) - odr = ODR(data, model, beta0=[1e-3]) - output = odr.run() - coeffsQ[i] = output.beta[0] - errorsQ[i] = output.sd_beta[0] - redchi2sQ[i] = output.sum_square / (len(x) - 2) - except: - coeffsQ[i] = np.nan - errorsQ[i] = np.nan - redchi2sQ[i] = np.nan - - #fit U - U_wnl = np.nanmean(np.abs(U_fs[i][wn_mask])) - x = np.real(I_fs[i])[fit_mask] - y = np.real(U_fs[i])[fit_mask] - sx = U_wnl / np.sqrt(2) * np.ones_like(x) - sy = U_wnl * np.ones_like(y) - try: - data = RealData(x=x, - y=y, - sx=sx, - sy=sy) - odr = ODR(data, model, beta0=[1e-3]) - output = odr.run() - coeffsU[i] = output.beta[0] - errorsU[i] = output.sd_beta[0] - redchi2sU[i] = output.sum_square / (len(x) - 2) - except: - coeffsU[i] = np.nan - errorsU[i] = np.nan - redchi2sU[i] = np.nan - - out_aman = core.AxisManager(aman.dets, aman.samps) - out_aman.wrap('coeffsQ', coeffsQ, [(0, 'dets')]) - out_aman.wrap('errorsQ', errorsQ, [(0, 'dets')]) - out_aman.wrap('redchi2sQ', redchi2sQ, [(0, 'dets')]) - out_aman.wrap('coeffsU', coeffsU, [(0, 'dets')]) - out_aman.wrap('errorsU', errorsU, [(0, 'dets')]) - out_aman.wrap('redchi2sU', redchi2sU, [(0, 'dets')]) + if ML_fit: + x = np.real(I_fs[:, fit_mask]) + yQ = np.real(Q_fs[:, fit_mask]) + yU = np.real(U_fs[:, fit_mask]) + coeffsQ = np.sum(x * yQ, axis=1) / np.sum(x**2, axis=1) + coeffsU = np.sum(x * yU, axis=1) / np.sum(x**2, axis=1) + # TODO: generalize this when using on data to apply T2P cuts. + errorsQ = np.zeros_like(coeffsQ) + errorsU = np.zeros_like(coeffsU) + redchi2sQ = np.zeros_like(coeffsQ) + redchi2sU = np.zeros_like(coeffsU) + else: + coeffsQ = np.zeros(aman.dets.count) + errorsQ = np.zeros(aman.dets.count) + redchi2sQ = np.zeros(aman.dets.count) + coeffsU = np.zeros(aman.dets.count) + errorsU = np.zeros(aman.dets.count) + redchi2sU = np.zeros(aman.dets.count) + + + wn_mask = (wn_freq_range[0] < freqs) & (freqs < wn_freq_range[1]) + + def leakage_model(B, x): + return B[0] * x + + model = Model(leakage_model) + + for i in range(aman.dets.count): + # fit Q + Q_wnl = np.nanmean(np.abs(Q_fs[i][wn_mask])) + x = np.real(I_fs[i])[fit_mask] + y = np.real(Q_fs[i])[fit_mask] + sx = Q_wnl / np.sqrt(2) * np.ones_like(x) + sy = Q_wnl * np.ones_like(y) + try: + data = RealData(x=x, + y=y, + sx=sx, + sy=sy) + odr = ODR(data, model, beta0=[1e-3]) + output = odr.run() + coeffsQ[i] = output.beta[0] + errorsQ[i] = output.sd_beta[0] + redchi2sQ[i] = output.sum_square / (len(x) - 2) + except: + coeffsQ[i] = np.nan + errorsQ[i] = np.nan + redchi2sQ[i] = np.nan + + #fit U + U_wnl = np.nanmean(np.abs(U_fs[i][wn_mask])) + x = np.real(I_fs[i])[fit_mask] + y = np.real(U_fs[i])[fit_mask] + sx = U_wnl / np.sqrt(2) * np.ones_like(x) + sy = U_wnl * np.ones_like(y) + try: + data = RealData(x=x, + y=y, + sx=sx, + sy=sy) + odr = ODR(data, model, beta0=[1e-3]) + output = odr.run() + coeffsU[i] = output.beta[0] + errorsU[i] = output.sd_beta[0] + redchi2sU[i] = output.sum_square / (len(x) - 2) + except: + coeffsU[i] = np.nan + errorsU[i] = np.nan + redchi2sU[i] = np.nan + + out_aman = core.AxisManager(aman.dets, aman.samps) + out_aman.wrap('coeffsQ', coeffsQ, [(0, 'dets')]) + out_aman.wrap('errorsQ', errorsQ, [(0, 'dets')]) + out_aman.wrap('redchi2sQ', redchi2sQ, [(0, 'dets')]) + out_aman.wrap('coeffsU', coeffsU, [(0, 'dets')]) + out_aman.wrap('errorsU', errorsU, [(0, 'dets')]) + out_aman.wrap('redchi2sU', redchi2sU, [(0, 'dets')]) if subtract_sig: subtract_t2p(aman, out_aman) From 20eaa20d565db33eabdec9ed800abe65070e611b Mon Sep 17 00:00:00 2001 From: Adrien La Posta Date: Tue, 24 Mar 2026 12:37:46 +0000 Subject: [PATCH 10/12] add Q std + red chi-sq --- sotodlib/tod_ops/t2pleakage.py | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/sotodlib/tod_ops/t2pleakage.py b/sotodlib/tod_ops/t2pleakage.py index 3e31d6914..d80f465e3 100644 --- a/sotodlib/tod_ops/t2pleakage.py +++ b/sotodlib/tod_ops/t2pleakage.py @@ -360,11 +360,15 @@ def get_t2p_coeffs_in_freq(aman, T_sig_name='dsT', Q_sig_name='demodQ', U_sig_na yU = np.real(U_fs[:, fit_mask]) coeffsQ = np.sum(x * yQ, axis=1) / np.sum(x**2, axis=1) coeffsU = np.sum(x * yU, axis=1) / np.sum(x**2, axis=1) - # TODO: generalize this when using on data to apply T2P cuts. + + stdQ = np.sqrt(np.sum((yQ - coeffsQ[:, np.newaxis] * x)**2, axis=1) / (x.shape[1] - 1)) + stdU = np.sqrt(np.sum((yU - coeffsU[:, np.newaxis] * x)**2, axis=1) / (x.shape[1] - 1)) + errorsQ = np.zeros_like(coeffsQ) errorsU = np.zeros_like(coeffsU) - redchi2sQ = np.zeros_like(coeffsQ) - redchi2sU = np.zeros_like(coeffsU) + + redchi2sQ = np.sum((yQ - coeffsQ[:, np.newaxis] * x) ** 2 / stdQ[:, np.newaxis] ** 2, axis=1) / (x.shape[1] - 1) + redchi2sU = np.sum((yU - coeffsU[:, np.newaxis] * x) ** 2 / stdU[:, np.newaxis] ** 2, axis=1) / (x.shape[1] - 1) else: coeffsQ = np.zeros(aman.dets.count) errorsQ = np.zeros(aman.dets.count) @@ -424,13 +428,13 @@ def leakage_model(B, x): errorsU[i] = np.nan redchi2sU[i] = np.nan - out_aman = core.AxisManager(aman.dets, aman.samps) - out_aman.wrap('coeffsQ', coeffsQ, [(0, 'dets')]) - out_aman.wrap('errorsQ', errorsQ, [(0, 'dets')]) - out_aman.wrap('redchi2sQ', redchi2sQ, [(0, 'dets')]) - out_aman.wrap('coeffsU', coeffsU, [(0, 'dets')]) - out_aman.wrap('errorsU', errorsU, [(0, 'dets')]) - out_aman.wrap('redchi2sU', redchi2sU, [(0, 'dets')]) + out_aman = core.AxisManager(aman.dets, aman.samps) + out_aman.wrap('coeffsQ', coeffsQ, [(0, 'dets')]) + out_aman.wrap('errorsQ', errorsQ, [(0, 'dets')]) + out_aman.wrap('redchi2sQ', redchi2sQ, [(0, 'dets')]) + out_aman.wrap('coeffsU', coeffsU, [(0, 'dets')]) + out_aman.wrap('errorsU', errorsU, [(0, 'dets')]) + out_aman.wrap('redchi2sU', redchi2sU, [(0, 'dets')]) if subtract_sig: subtract_t2p(aman, out_aman) From b8a241b4d7e20b301c1dad8da170bebbbf8b5b59 Mon Sep 17 00:00:00 2001 From: Adrien La Posta Date: Tue, 24 Mar 2026 12:39:10 +0000 Subject: [PATCH 11/12] update docstring --- sotodlib/tod_ops/t2pleakage.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/sotodlib/tod_ops/t2pleakage.py b/sotodlib/tod_ops/t2pleakage.py index d80f465e3..aa2642239 100644 --- a/sotodlib/tod_ops/t2pleakage.py +++ b/sotodlib/tod_ops/t2pleakage.py @@ -337,6 +337,10 @@ def get_t2p_coeffs_in_freq(aman, T_sig_name='dsT', Q_sig_name='demodQ', U_sig_na ML_fit : bool Whether to use the ML solution (analytic) to fit for coefficients. Default is False, which uses ODR to perform a linear fit. + We are assuming a linear model: y_i = b * x_i + eps_i, where eps_i is white noise. + The log likelihood for this model is: log L(b) = - (1 / (2*sigma^2)) * sum_i (y_i - b*x_i)^2 + const + For which the maximum likelihood estimator for b is: b_hat = sum_i (x_i * y_i) / sum_i (x_i^2) + This is equivalent to Ordinary Least Square (OLS) regression with no intercept. Returns ------- From cace6044269b496b629c87a7600b542243a19f5c Mon Sep 17 00:00:00 2001 From: Adrien La Posta Date: Tue, 24 Mar 2026 15:08:41 +0000 Subject: [PATCH 12/12] ML errors for sims --- sotodlib/tod_ops/t2pleakage.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sotodlib/tod_ops/t2pleakage.py b/sotodlib/tod_ops/t2pleakage.py index aa2642239..972105b29 100644 --- a/sotodlib/tod_ops/t2pleakage.py +++ b/sotodlib/tod_ops/t2pleakage.py @@ -368,8 +368,8 @@ def get_t2p_coeffs_in_freq(aman, T_sig_name='dsT', Q_sig_name='demodQ', U_sig_na stdQ = np.sqrt(np.sum((yQ - coeffsQ[:, np.newaxis] * x)**2, axis=1) / (x.shape[1] - 1)) stdU = np.sqrt(np.sum((yU - coeffsU[:, np.newaxis] * x)**2, axis=1) / (x.shape[1] - 1)) - errorsQ = np.zeros_like(coeffsQ) - errorsU = np.zeros_like(coeffsU) + errorsQ = stdQ / np.sqrt(np.sum(x**2, axis=1)) + errorsU = stdU / np.sqrt(np.sum(x**2, axis=1)) redchi2sQ = np.sum((yQ - coeffsQ[:, np.newaxis] * x) ** 2 / stdQ[:, np.newaxis] ** 2, axis=1) / (x.shape[1] - 1) redchi2sU = np.sum((yU - coeffsU[:, np.newaxis] * x) ** 2 / stdU[:, np.newaxis] ** 2, axis=1) / (x.shape[1] - 1)