diff --git a/sotodlib/coords/det_match.py b/sotodlib/coords/det_match.py index 97f7596f4..e13410ccf 100644 --- a/sotodlib/coords/det_match.py +++ b/sotodlib/coords/det_match.py @@ -613,14 +613,14 @@ class MatchParams: unassigned_slots: int = 1000 freq_offset_mhz: float = 0.0 freq_width: float = 2. - dist_width: float =0.01 + dist_width: float = 0.01 unmatched_good_res_pen: float = 10. good_res_qi_thresh: float = 100e3 enforce_pointing_reqs: bool = False assigned_bg_unmatched_pen: float = 100000 unassigned_bg_unmatched_pen: float = 10000 - assigned_bg_mismatch_pen: float = 100000 + assigned_bg_mismatch_pen: float = np.inf #100000 unassigned_bg_mismatch_pen: float = 1 @@ -696,9 +696,11 @@ def __init__(self, src: ResSet, dst: ResSet, match_pars: Optional[MatchParams]=N self.matching, self.merged = self._match() self.stats = self.get_stats() + # We're matching dst onto src so design onto pointing or + # pointing onto tune def _get_biadjacency_matrix(self) -> np.ndarray: - src_arr = self.src.as_array() - dst_arr = self.dst.as_array() + src_arr = self.src.as_array() # pointing or tuneset + dst_arr = self.dst.as_array() # design or pointing mat = np.zeros((len(self.src), len(self.dst)), dtype=float) @@ -711,14 +713,17 @@ def _get_biadjacency_matrix(self) -> np.ndarray: # pointing data. For these types of detectors, the cost is left # untouched. + # Has pointing (design and pointing) src_has_pointing = np.isfinite(src_arr['xi']) & np.isfinite(src_arr['eta']) dst_has_pointing = np.isfinite(dst_arr['xi']) & np.isfinite(dst_arr['eta']) + # Do not match NC src_no_match = np.isin(src_arr['det_type'], ['NC']) dst_no_match = np.isin(dst_arr['det_type'], ['NC']) mat[src_no_match, :] = np.inf mat[:, dst_no_match] = np.inf + # Don't want pointing matches for these detector types src_pointing_forbidden = np.isin(src_arr['det_type'], ['UNRT', 'SQID', 'BARE']) dst_pointing_forbidden = np.isin(dst_arr['det_type'], ['UNRT', 'SQID', 'BARE']) m = src_pointing_forbidden[:, None] & dst_has_pointing[None, :] @@ -726,6 +731,7 @@ def _get_biadjacency_matrix(self) -> np.ndarray: m = src_has_pointing[:, None] & dst_pointing_forbidden[None, :] mat[m] = np.inf + # OPTC detectors should have pointing src_pointing_required = np.isin(src_arr['det_type'], ['OPTC']) dst_pointing_required = np.isin(dst_arr['det_type'], ['OPTC']) m = src_pointing_required[:, None] & (~dst_has_pointing[None, :]) @@ -733,6 +739,16 @@ def _get_biadjacency_matrix(self) -> np.ndarray: m = (~src_has_pointing[:, None]) & dst_pointing_required[None, :] mat[m] = np.inf + # These should always have BG = -1 + src_unassigned_type = np.isin(src_arr['det_type'], ['UNRT', 'SQID', 'BARE']) + dst_unassigned_type = np.isin(dst_arr['det_type'], ['UNRT', 'SQID', 'BARE']) + # pointing or tuneset have unassigned type (shouldn't ever happen) + m = src_unassigned_type[:, None] & (dst_arr['bg'][None, :] !=-1) + mat[m] = np.inf + # Design or pointing have unassigned type (should happen) + m = dst_unassigned_type[None, :] & (src_arr['bg'][:, None] !=-1) + mat[m] = np.inf + # Frequency offset df = src_arr['res_freq'][:, None] - dst_arr['res_freq'][None, :] df -= self.match_pars.freq_offset_mhz @@ -740,10 +756,12 @@ def _get_biadjacency_matrix(self) -> np.ndarray: # BG mismatch bgs_mismatch = src_arr['bg'][:, None] != dst_arr['bg'][None, :] - bgs_unassigned = (src_arr['bg'][:, None] == 1) | (dst_arr['bg'][None, :] == -1) - + # If either BG is unassigned + bgs_unassigned = (src_arr['bg'][:, None] == -1) | (dst_arr['bg'][None, :] == -1) + # If the BG is mismatched and one of them is unassigned m = bgs_mismatch & bgs_unassigned mat[m] += self.match_pars.unassigned_bg_mismatch_pen + # If the BG is mismatched and none of them are unassigned (shouldn't happen) m = bgs_mismatch & (~bgs_unassigned) mat[m] += self.match_pars.assigned_bg_mismatch_pen @@ -759,6 +777,7 @@ def _get_biadjacency_matrix(self) -> np.ndarray: return mat + def _get_unassigned_costs(self, rs): ra = rs.as_array() @@ -779,6 +798,7 @@ def _get_unassigned_costs(self, rs): def _match(self): nside = max(len(self.src), len(self.dst)) + self.match_pars.unassigned_slots + # Keep this square so all resonators are included in final matching mat_full = np.zeros((nside, nside), dtype=float) with warnings.catch_warnings(): @@ -798,14 +818,28 @@ def _match(self): self.matching = np.array(linear_sum_assignment(mat_full)) self.matching_cost = mat_full[self.matching[0], self.matching[1]].sum() + # ensure matching and match_idx are reset for r1, r2 in self.get_match_iter(include_unmatched=True): if r1 is None: r2.matched = 0 + r2.match_idx = -1 continue if r2 is None: r1.matched = 0 + r1.match_idx = -1 continue + r1.matched = 0 + r1.match_idx = -1 + r2.matched = 0 + r2.match_idx = -1 + for r1, r2 in self.get_match_iter(include_unmatched=True): + if r1 is None: + r2.matched = 0 + continue + if r2 is None: + r1.matched = 0 + continue r1.matched = 1 r1.match_idx = r2.idx r2.matched = 1 diff --git a/sotodlib/site_pipeline/update_det_match.py b/sotodlib/site_pipeline/update_det_match.py index 08fca4ef9..dfd6f3196 100644 --- a/sotodlib/site_pipeline/update_det_match.py +++ b/sotodlib/site_pipeline/update_det_match.py @@ -80,7 +80,7 @@ class UpdateDetMatchesConfig: Time in seconds before a failed detset will be added to the cache. This is to prevent new detsets still acquiring data from being added right away. - + Attributes ------------- freq_offsets : Optional[np.ndarray] @@ -109,7 +109,7 @@ def __post_init__(self): if self.wafer_map_path is None: self.wafer_map_path = os.path.join( self.site_pipeline_root, 'shared/detmapping/wafer_map.yaml') - + if self.freq_offset_range_args is not None: self.freq_offsets = np.arange(*self.freq_offset_range_args) else: @@ -120,12 +120,12 @@ def __post_init__(self): if not os.path.exists(self.results_path): raise FileNotFoundError(f"Results dir does not exist: {self.results_path}") - + allowed_solution_types = ['kaiwen_handmade', 'resonator_set'] if self.solution_type not in allowed_solution_types: raise ValueError( f"Solution type ({self.solution_type}) must be a member of: {allowed_solution_types}") - + if self.solution_type == 'resonator_set': if self.resonator_set_dir is None: raise ValueError("Must specify resonator_set_dir for solution_type='resonator_set'") @@ -139,7 +139,7 @@ def __init__(self, cfg: UpdateDetMatchesConfig): self.detcal_db = None with open(self.cfg.wafer_map_path, 'r') as f: self.wafer_map = yaml.safe_load(f) - + self.ufm_to_fp_file = os.path.join( cfg.site_pipeline_root, 'shared/focalplane/ufm_to_fp.yaml') @@ -161,7 +161,7 @@ def __init__(self, cfg: UpdateDetMatchesConfig): if self.detcal_db is None: raise Exception( f"Could not find detcal metadata entry with name: {cfg.detcal_meta_name}") - + self.failed_detset_cache_path = os.path.join(cfg.results_path, 'failed_detsets.yaml') self.match_dir = os.path.join(cfg.results_path, 'matches') if not os.path.exists(self.match_dir): @@ -212,7 +212,7 @@ def add_to_failed_cache(cache_file, detset, msg, cfg: UpdateDetMatchesConfig): with open(cache_file, 'r') as f: x = yaml.safe_load(f) else: - x = {} + x = {} x[str(detset)] = str(msg) with open(cache_file, 'w') as f: @@ -244,7 +244,7 @@ def run_match_aman(runner: Runner, aman, detset, wafer_slot=None): match_pars=match_pars, ) match_pars.freq_offset_mhz = opt_freq - match = det_match.Match(rs0, rs1, match_pars=match_pars, + match = det_match.Match(rs0, rs1, match_pars=match_pars, apply_dst_pointing=runner.cfg.apply_solution_pointing) return match @@ -267,7 +267,7 @@ def run_match(runner: Runner, detset: str) -> bool: obsids_with_cal = set(runner.detcal_db.get_entries(['dataset'])['dataset']) obs_dset = {r[0] for r in cur} obs_ids = sorted(list( - obs_all.intersection(obs_dset).intersection(obsids_with_cal)), + obs_all.intersection(obs_dset).intersection(obsids_with_cal)), key=lambda s:s.split('_')[1])[::-1] if len(obs_ids) == 0: add_to_failed_cache( @@ -280,8 +280,11 @@ def run_match(runner: Runner, detset: str) -> bool: book_dir = os.path.dirname(runner.ctx.obsfiledb.get_files(obs_id)[detset][0][0]) book_idx_file = os.path.join(book_dir, 'M_index.yaml') - with open(book_idx_file, 'r') as f: - book_idx = yaml.safe_load(f) + try: + with open(book_idx_file, 'r') as f: + book_idx = yaml.safe_load(f) + except Exception as e: + return False aman = runner.ctx.get_meta(obs_id, ignore_missing=True) finished_detsets = set([os.path.splitext(f)[0] for f in os.listdir(runner.match_dir)]) @@ -299,7 +302,7 @@ def run_match(runner: Runner, detset: str) -> bool: stream_id = aman.det_info.stream_id[aman.det_info.detset == ds][0] # Try to get wafer slot info from book idx if 'wafer_slots' in book_idx: - for ws in book_idx['wafer_slots']: + for ws in book_idx['wafer_slots']: if ws['stream_id'] == stream_id: wafer_slot = ws['wafer_slot'] break @@ -311,11 +314,14 @@ def run_match(runner: Runner, detset: str) -> bool: raise Exception("Could not find wafer-slot") else: wafer_slot = None - - match = run_match_aman(runner, aman, ds, wafer_slot=wafer_slot) - fpath = os.path.join(runner.match_dir, f"{ds}.h5") - match.save(fpath) - logger.info(f"Saved match to file: {fpath}") + try: + match = run_match_aman(runner, aman, ds, wafer_slot=wafer_slot) + fpath = os.path.join(runner.match_dir, f"{ds}.h5") + match.save(fpath) + logger.info(f"Saved match to file: {fpath}") + except Exception as e: + logger.error(f"{ds} failed with {e}") + continue return True @@ -338,8 +344,7 @@ def scan_for_freq_offset(rs0, rs1, freq_offsets, match_pars=None, show_pb=True): else: match_pars = deepcopy(match_pars) freq_offsets = np.array(freq_offsets) - - rs0 = deepcopy(rs0) + rs0 = deepcopy(rs0) rs1 = deepcopy(rs1) costs = np.full_like(freq_offsets, np.nan) @@ -389,7 +394,7 @@ def add_to_db(arr, db_path, h5_path, detset, write_relpath=True): }, h5_path) else: logger.warning(f"Dataset {detset} already exists in db: {db_path}") - + write_relpath = runner.cfg.write_relpath add_to_db(ra, det_match_idx, det_match_h5, detset, write_relpath=write_relpath)