Skip to content

Det match updates #1515

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.

Sign up for GitHub

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

Closed
mmccrackan wants to merge 1 commit into from
Closed

Det match updates #1515

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
46 changes: 40 additions & 6 deletions sotodlib/coords/det_match.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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


Expand Down Expand Up @@ -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)

Expand All @@ -711,39 +713,55 @@ 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, :]
mat[m] = np.inf
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, :])
mat[m] = np.inf
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
mat += np.exp((np.abs(df / self.match_pars.freq_width)) ** 2)

# 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

Expand All @@ -759,6 +777,7 @@ def _get_biadjacency_matrix(self) -> np.ndarray:

return mat


def _get_unassigned_costs(self, rs):
ra = rs.as_array()

Expand All @@ -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():
Expand All @@ -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
Expand Down
45 changes: 25 additions & 20 deletions sotodlib/site_pipeline/update_det_match.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down Expand Up @@ -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:
Expand All @@ -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'")
Expand All @@ -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')

Expand All @@ -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):
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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

Expand All @@ -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(
Expand All @@ -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)])
Expand All @@ -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
Expand All @@ -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

Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand Down