From e2f2eedb3698474f568d7b4d9d5a5c139be26c28 Mon Sep 17 00:00:00 2001 From: Javier Gonzalez Date: Tue, 17 Mar 2026 14:09:43 -0400 Subject: [PATCH 01/14] update pre-commit --- .pre-commit-config.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index f82d48f..7343548 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,7 +1,7 @@ repos: - repo: https://github.com/astral-sh/ruff-pre-commit # Ruff version. - rev: v0.5.1 + rev: v0.14.13 hooks: # Run the linter. - id: ruff From 2bc35fe73fe00efd698f4641cfb3e0066bf5c88f Mon Sep 17 00:00:00 2001 From: Javier Gonzalez Date: Tue, 17 Mar 2026 14:11:57 -0400 Subject: [PATCH 02/14] separate the obs-stats part from the star-specific part in update_mag_supplement.get_agasc_id_stats --- agasc/supplement/magnitudes/mag_estimate.py | 523 ++++++++++++-------- 1 file changed, 313 insertions(+), 210 deletions(-) diff --git a/agasc/supplement/magnitudes/mag_estimate.py b/agasc/supplement/magnitudes/mag_estimate.py index 6f623ea..23f41db 100644 --- a/agasc/supplement/magnitudes/mag_estimate.py +++ b/agasc/supplement/magnitudes/mag_estimate.py @@ -481,7 +481,9 @@ def get_telemetry(obs): return telem -def get_telemetry_by_agasc_id(agasc_id, obsid=None, ignore_exceptions=False): +def get_telemetry_by_agasc_id( + agasc_id, obsid=None, ignore_exceptions=False, as_table=True +): """ Get all telemetry relevant for the magnitude estimation, given an AGASC ID. @@ -491,10 +493,13 @@ def get_telemetry_by_agasc_id(agasc_id, obsid=None, ignore_exceptions=False): :param obsid: int (optional) :param ignore_exceptions: bool if True, any exception is ignored. Useful in some cases. Default is False. + :param as_table: bool + if True, the telemetry is returned as an astropy Table. If False, it is returned as a list + of dicts, one per observation. Default is True. :return: dict """ logger.debug(f" Getting telemetry for AGASC ID={agasc_id}") - star_obs_catalogs.load() + star_obs_catalogs.load() # should this be here? if obsid is None: obs = star_obs_catalogs.STARS_OBS[ (star_obs_catalogs.STARS_OBS["agasc_id"] == agasc_id) @@ -505,25 +510,48 @@ def get_telemetry_by_agasc_id(agasc_id, obsid=None, ignore_exceptions=False): & (star_obs_catalogs.STARS_OBS["obsid"] == obsid) ] obs.sort("mp_starcat_time") + return get_telemetry_by_observations( + obs, ignore_exceptions=ignore_exceptions, as_table=as_table + ) + +def get_telemetry_by_observations(observations, ignore_exceptions=False, as_table=True): telem = [] - for _i, o in enumerate(obs): + for obs in observations: + agasc_id = obs["agasc_id"] try: - t = Table(get_telemetry(o)) - t["obsid"] = o["obsid"] - t["agasc_id"] = agasc_id + t = get_telemetry(obs) + t["obsid"] = obs["obsid"] * np.ones(len(t["times"]), dtype=int) + t["agasc_id"] = agasc_id * np.ones(len(t["times"]), dtype=int) telem.append(t) + except MagStatsException as exc: + if ignore_exceptions: + telem.append(dict(exc)) + else: + logger.info(f"{agasc_id=}, obsid={obs['obsid']} failed") + logger.info(f"{exc.exception['name']} {exc.exception['value']}") + for step in exc.exception["traceback"]: + logger.info(step) except Exception: - if not ignore_exceptions: - logger.info(f"{agasc_id=}, obsid={o['obsid']} failed") - exc_type, exc_value, exc_traceback = sys.exc_info() - trace = traceback.extract_tb(exc_traceback) - logger.info(f"{exc_type.__name__} {exc_value}") - for step in trace: - logger.info(f" in {step.filename}:{step.lineno}/{step.name}:") - logger.info(f" {step.line}") + exc = MagStatsException( + msg="Unknown", + agasc_id=obs["agasc_id"], + obsid=obs["obsid"], + mp_starcat_time=obs["mp_starcat_time"], + ) + if ignore_exceptions: + telem.append(dict(exc)) + else: + logger.info(f"{agasc_id=}, obsid={obs['obsid']} failed") + logger.info(f"{exc.exception['name']} {exc.exception['value']}") + for step in exc.exception["traceback"]: + logger.info(step) raise - return vstack(telem) if telem else [] + + if telem: + return vstack([Table(tel) for tel in telem]) if as_table else telem + + return [] def add_obs_info(telem, obs_stats): @@ -797,6 +825,7 @@ def get_obs_stats(obs, telem=None): "mag_aca_err": star["MAG_ACA_ERR"] / 100, "row": obs["row"], "col": obs["col"], + "obs_ok": False, } ) @@ -862,6 +891,13 @@ def get_obs_stats(obs, telem=None): f"f_mag_est_ok={stats['f_mag_est_ok']:.3f}, f_dr3={stats['f_dr3']:.3f}, " f"mag={stats['mag_obs']:.2f}" ) + + stats["obs_ok"] = ( + (stats["n"] > 10) + & (stats["f_mag_est_ok"] > 0.3) + & (stats["lf_variability_100s"] < 1) + ) + return stats @@ -1168,6 +1204,55 @@ def get_agasc_id_stats(agasc_id, obs_status_override=None, tstop=None): ] star_obs.sort("mp_starcat_time") + n_obsids = len(star_obs) + + all_telem = get_telemetry_by_observations( + star_obs, ignore_exceptions=True, as_table=False + ) + stats, failures = get_multi_obs_stats( + star_obs, obs_status_override=obs_status_override, telem=all_telem + ) + + # combine magnitude estimates using a weighted mean + weighted_mean = get_weighted_mean(stats) + stats["w"] = weighted_mean["weights"] + stats["mean_corrected"] = weighted_mean["mean_corrected"] + stats["weighted_mean"] = weighted_mean["weighted_mean"] + + star = get_star(agasc_id, use_supplement=False) + + # still need to check that this is the same as before + last_obs_time = CxoTime(stats["mp_starcat_time"][-1]).cxcsec + + logger.debug(" identifying outlying observations...") + for s, t in zip(stats, all_telem, strict=True): + if s["no_telem"] or s["excluded"]: + continue + t["obs_ok"] = np.ones_like(t["mag_est_ok"], dtype=bool) * s["obs_ok"] + logger.debug( + " identifying outlying observations " + f"(OBSID={s['obsid']}, mp_starcat_time={s['mp_starcat_time']})" + ) + t["obs_outlier"] = np.zeros_like(t["mag_est_ok"]) + if np.any(t["mag_est_ok"]) and s["f_mag_est_ok"] > 0 and s["obs_ok"]: + iqr = s["q75"] - s["q25"] + t["obs_outlier"] = ( + t["mag_est_ok"] + & (iqr > 0) + & ( + (t["mags"] < s["q25"] - 1.5 * iqr) + | (t["mags"] > s["q75"] + 1.5 * iqr) + ) + ) + + all_telem = [ + Table(t) + for s, t in zip(stats, all_telem, strict=True) + if not s["excluded"] and not s["no_telem"] + ] + if len(all_telem) > 0: + all_telem = vstack(all_telem) + # this is the default result, if nothing gets calculated result = { "last_obs_time": 0, @@ -1196,8 +1281,8 @@ def get_agasc_id_stats(agasc_id, obs_status_override=None, tstop=None): "sigma_plus": 0, "mean": 0, "std": 0, - "mag_weighted_mean": 0, - "mag_weighted_std": 0, + "mag_weighted_mean": weighted_mean["mag_weighted_mean"], + "mag_weighted_std": weighted_mean["mag_weighted_std"], "t_mean": 0, "t_std": 0, "n_outlier": 0, @@ -1212,13 +1297,31 @@ def get_agasc_id_stats(agasc_id, obs_status_override=None, tstop=None): "selected_rtol": False, "selected_mag_aca_err": False, "selected_color": False, - "f_mag_est_ok": 0, - "f_mag_est_ok_3": 0, - "f_mag_est_ok_5": 0, - "f_ok_3": 0, - "f_ok_5": 0, + "f_mag_est_ok": 0.0, + "f_mag_est_ok_3": 0.0, + "f_mag_est_ok_5": 0.0, + "f_ok_3": 0.0, + "f_ok_5": 0.0, } + # this can be moved up + result.update( + { + "color": star["COLOR1"], + "last_obs_time": last_obs_time, + "mag_aca": star["MAG_ACA"], + "mag_aca_err": star["MAG_ACA_ERR"] / 100, + "mag_obs_err": min_mag_obs_err, + "n_obsids_fail": len(failures), + "n": len(all_telem), + "n_obsids_suspect": np.count_nonzero(stats["obs_suspect"]), + "n_obsids": n_obsids, + "n_obsids_ok": np.count_nonzero(stats["obs_ok"]), + "n_no_mag": np.count_nonzero((~stats["obs_ok"])) + + np.count_nonzero(stats["f_mag_est_ok"][stats["obs_ok"]] < 0.3), + } + ) + for tag in ["dr3", "dbox5"]: result.update( { @@ -1237,127 +1340,8 @@ def get_agasc_id_stats(agasc_id, obs_status_override=None, tstop=None): } ) - n_obsids = len(star_obs) - - # exclude star_obs that are in obs_status_override with status != 0 - excluded_obs = np.array( - [ - ( - (oi, ai) in obs_status_override - and obs_status_override[(oi, ai)]["status"] != 0 - ) - for oi, ai in star_obs[["mp_starcat_time", "agasc_id"]] - ] - ) - if np.any(excluded_obs): - logger.debug( - " Excluding observations flagged in obs-status table: " - f"{list(star_obs[excluded_obs]['obsid'])}" - ) - - included_obs = np.array( - [ - ( - (oi, ai) in obs_status_override - and obs_status_override[(oi, ai)]["status"] == 0 - ) - for oi, ai in star_obs[["mp_starcat_time", "agasc_id"]] - ] - ) - if np.any(included_obs): - logger.debug( - " Including observations marked OK in obs-status table: " - f"{list(star_obs[included_obs]['obsid'])}" - ) - - failures = [] - all_telem = [] - stats = [] - last_obs_time = 0 - for i, obs in enumerate(star_obs): - oi, ai = obs["mp_starcat_time", "agasc_id"] - comment = "" - if (oi, ai) in obs_status_override: - status = obs_status_override[(oi, ai)] - logger.debug( - f" overriding status for (AGASC ID {ai}, starcat time {oi}): " - f"{status['status']}, {status['comments']}" - ) - comment = status["comments"] - try: - last_obs_time = CxoTime(obs["mp_starcat_time"]).cxcsec - telem = Table(get_telemetry(obs)) - obs_stat = get_obs_stats(obs, telem={k: telem[k] for k in telem.colnames}) - obs_stat.update( - { - "obs_ok": included_obs[i] - | ( - ~excluded_obs[i] - & (obs_stat["n"] > 10) - & (obs_stat["f_mag_est_ok"] > 0.3) - & (obs_stat["lf_variability_100s"] < 1) - ), - "obs_suspect": False, - "obs_fail": False, - "comments": comment, - } - ) - all_telem.append(telem) - stats.append(obs_stat) - - if not obs_stat["obs_ok"] and not excluded_obs[i]: - obs_stat["obs_suspect"] = True - failures.append( - dict( - MagStatsException( - msg="Suspect observation", - agasc_id=obs["agasc_id"], - obsid=obs["obsid"], - mp_starcat_time=obs["mp_starcat_time"], - ) - ) - ) - except MagStatsException as e: - # this except branch deals with exceptions thrown by get_telemetry - all_telem.append(None) - # length-zero telemetry short-circuits any new call to get_telemetry - obs_stat = get_obs_stats(obs, telem=[]) - obs_stat.update( - { - "obs_ok": False, - "obs_suspect": False, - "obs_fail": e.failed, - "comments": comment if excluded_obs[i] else f"Error: {e.msg}.", - } - ) - stats.append(obs_stat) - if e.failed and not excluded_obs[i]: - logger.debug( - f" Error in get_agasc_id_stats({agasc_id=}," - f" obsid={obs['obsid']}): {e}" - ) - failures.append(dict(e)) - - stats = Table(stats) - stats["w"] = np.nan - stats["mean_corrected"] = np.nan - stats["weighted_mean"] = np.nan - - star = get_star(agasc_id, use_supplement=False) - - result.update( - { - "last_obs_time": last_obs_time, - "mag_aca": star["MAG_ACA"], - "mag_aca_err": star["MAG_ACA_ERR"] / 100, - "color": star["COLOR1"], - "n_obsids_fail": len(failures), - "n_obsids_suspect": np.count_nonzero(stats["obs_suspect"]), - "n_obsids": n_obsids, - } - ) - - if not np.any(~excluded_obs): + # Check requirements for calculating metrics. + if np.all(stats["excluded"]): # this is a new field # this happens when all observations have been flagged as not OK a priory (obs-status). logger.debug( f" Skipping star in get_agasc_id_stats({agasc_id=})." @@ -1365,37 +1349,16 @@ def get_agasc_id_stats(agasc_id, obs_status_override=None, tstop=None): ) return result, stats, failures - # add failed observations to the list of excluded observations - excluded_obs += np.array([t is None for t in all_telem]) - # and remove all excluded observations from all_telem - all_telem = [t for i, t in enumerate(all_telem) if not excluded_obs[i]] - - if len(all_telem) == 0: - # and we reach here if some observations were not flagged as bad, but all failed. + if np.count_nonzero(stats["excluded"] | stats["no_telem"]) == len(stats): + # and we reach here if some observations were not flagged as bad, but + # there was some error getting telemetry. logger.debug(f" get_agasc_id_stats({agasc_id=}): There is no OK observation.") return result, stats, failures - logger.debug(" identifying outlying observations...") - for s, t in zip(stats[~excluded_obs], all_telem, strict=True): - t["obs_ok"] = np.ones_like(t["mag_est_ok"], dtype=bool) * s["obs_ok"] - logger.debug( - " identifying outlying observations " - f"(OBSID={s['obsid']}, mp_starcat_time={s['mp_starcat_time']})" - ) - t["obs_outlier"] = np.zeros_like(t["mag_est_ok"]) - if np.any(t["mag_est_ok"]) and s["f_mag_est_ok"] > 0 and s["obs_ok"]: - iqr = s["q75"] - s["q25"] - t["obs_outlier"] = ( - t["mag_est_ok"] - & (iqr > 0) - & ( - (t["mags"] < s["q25"] - 1.5 * iqr) - | (t["mags"] > s["q75"] + 1.5 * iqr) - ) - ) - all_telem = vstack([Table(t) for t in all_telem]) - n_total = len(all_telem) + if np.count_nonzero(all_telem["mag_est_ok"] & all_telem["obs_ok"]) < 10: + return result, stats, failures + # First calculate metrics that do not depend on centroid filters kalman = (all_telem["AOACASEQ"] == "KALM") & (all_telem["AOPCADMD"] == "NPNT") all_telem = all_telem[kalman] # non-npm/non-kalman are excluded n_kalman = len(all_telem) @@ -1406,23 +1369,6 @@ def get_agasc_id_stats(agasc_id, obs_status_override=None, tstop=None): sat_pix = (all_telem["AOACISP"] == "OK") & all_telem["obs_ok"] ion_rad = (all_telem["AOACIIR"] == "OK") & all_telem["obs_ok"] - f_mag_est_ok = np.count_nonzero(mag_est_ok) / len(mag_est_ok) - - result.update( - { - "mag_obs_err": min_mag_obs_err, - "n_obsids_ok": np.count_nonzero(stats["obs_ok"]), - "n_no_mag": np.count_nonzero((~stats["obs_ok"])) - + np.count_nonzero(stats["f_mag_est_ok"][stats["obs_ok"]] < 0.3), - "n": n_total, - "n_mag_est_ok": np.count_nonzero(mag_est_ok), - "f_mag_est_ok": f_mag_est_ok, - } - ) - - if result["n_mag_est_ok"] < 10: - return result, stats, failures - sigma_minus, q25, median, q75, sigma_plus = np.quantile( mags[mag_est_ok], [0.158, 0.25, 0.5, 0.75, 0.842] ) @@ -1431,37 +1377,16 @@ def get_agasc_id_stats(agasc_id, obs_status_override=None, tstop=None): outlier_2 = mag_est_ok & ((mags > q75 + 3 * iqr) | (mags < q25 - 3 * iqr)) outlier = all_telem["obs_outlier"] - # combine measurements using a weighted mean - obs_ok = stats["obs_ok"] - min_std = max(0.1, stats[obs_ok]["std"].min()) - stats["w"][obs_ok] = np.where( - stats["std"][obs_ok] != 0, 1.0 / stats["std"][obs_ok], 1.0 / min_std - ) - stats["mean_corrected"][obs_ok] = ( - stats["t_mean"][obs_ok] + stats["mag_correction"][obs_ok] - ) - stats["weighted_mean"][obs_ok] = ( - stats["mean_corrected"][obs_ok] * stats["w"][obs_ok] - ) - - mag_weighted_mean = stats[obs_ok]["weighted_mean"].sum() / stats[obs_ok]["w"].sum() - mag_weighted_std = np.sqrt( - ((stats[obs_ok]["mean"] - mag_weighted_mean) ** 2 * stats[obs_ok]["w"]).sum() - / stats[obs_ok]["w"].sum() - ) - result.update( { "agasc_id": agasc_id, "n_mag_est_ok": np.count_nonzero(mag_est_ok), - "f_mag_est_ok": f_mag_est_ok, + "f_mag_est_ok": np.count_nonzero(mag_est_ok) / len(mag_est_ok), "median": median, "sigma_minus": sigma_minus, "sigma_plus": sigma_plus, "mean": np.mean(mags[mag_est_ok]), "std": np.std(mags[mag_est_ok]), - "mag_weighted_mean": mag_weighted_mean, - "mag_weighted_std": mag_weighted_std, "t_mean": np.mean(mags[mag_est_ok & (~outlier)]), "t_std": np.std(mags[mag_est_ok & (~outlier)]), "n_outlier": np.count_nonzero(mag_est_ok & outlier), @@ -1474,8 +1399,9 @@ def get_agasc_id_stats(agasc_id, obs_status_override=None, tstop=None): } ) + # Calculate metrics that depend on centroid filters residual_ok = { - 3: all_telem["dr"] < 3, + 3: (np.sqrt(all_telem["dy"] ** 2 + all_telem["dz"] ** 2) < 3), 5: (np.abs(all_telem["dy"]) < 5) & (np.abs(all_telem["dz"]) < 5), } dr_tag = {3: "dr3", 5: "dbox5"} @@ -1521,6 +1447,7 @@ def get_agasc_id_stats(agasc_id, obs_status_override=None, tstop=None): } ) + # main results, which correspond to the "dbox5" residual cut result.update( { "mag_obs": result["t_mean_dbox5"], @@ -1547,3 +1474,179 @@ def get_agasc_id_stats(agasc_id, obs_status_override=None, tstop=None): logger.debug(f" stats for AGASC ID {agasc_id}: {stats['mag_obs'][0]}") return result, stats, failures + + +def get_weighted_mean(stats): + # combine measurements using a weighted mean + weights = np.nan * np.ones(len(stats)) + mean_corrected = np.nan * np.ones(len(stats)) + weighted_mean = np.nan * np.ones(len(stats)) + mag_weighted_mean = 0.0 + mag_weighted_std = 0.0 + + obs_ok = stats["obs_ok"] + if np.any(obs_ok): + min_std = max(0.1, stats[obs_ok]["std"].min()) + weights[obs_ok] = np.where( + stats["std"][obs_ok] != 0, 1.0 / stats["std"][obs_ok], 1.0 / min_std + ) + mean_corrected[obs_ok] = ( + stats["t_mean"][obs_ok] + stats["mag_correction"][obs_ok] + ) + weighted_mean[obs_ok] = mean_corrected[obs_ok] * weights[obs_ok] + + mag_weighted_mean = weighted_mean[obs_ok].sum() / weights[obs_ok].sum() + mag_weighted_std = np.sqrt( + ((stats[obs_ok]["mean"] - mag_weighted_mean) ** 2 * weights[obs_ok]).sum() + / weights[obs_ok].sum() + ) + + result = { + "weights": weights, + "mean_corrected": mean_corrected, + "weighted_mean": weighted_mean, + "mag_weighted_mean": mag_weighted_mean, + "mag_weighted_std": mag_weighted_std, + } + return result + + +def get_multi_obs_stats(star_obs, telem=None, obs_status_override=None): + """ + Get summary magnitude statistics for an AGASC ID. + + This function deals with several errors that can occur, and assigns various success/error flags + to each observation: + + - obs_ok: OK observations are considered OK based on criteria listed below, or are + have `status == 0` in `obs_status_override`. + - obs_suspect: observations are considered suspect if they have `status > 1` in + `obs_status_override`, and do not fulfill the criteria for being OK. + - obs_fail: observations where there was an exception getting the telemetry or calculating + the stats. + - no_telem: observations where there was an exception getting the telemetry. + + Criteria for an OK observation: + + - n > 10, where n is the total number of telemetry entries in the observation. + - f_mag_est_ok > 0.3, where f_mag_est_ok is the fraction of telemetry entries with + (AOACASEQ==KALM & AOACIIR==OK & AOPCADMD==NPNT & AOACFCT==TRAK). + - lf_variability_100s < 1, where lf_variability_100s is the largest change in the magnitude + smoothed using a 100s rolling window. + + :param star_obs: Table + Table of star observations. + :param telem: list + List of telemetry tables for each observation. Must be in the same order as star_obs. + This should generally be the result of get_telemetry_by_observations. + :param obs_status_override: dict. + Dictionary overriding the OK flag for specific observations. + Keys are (OBSID, AGASC ID) pairs, values are dictionaries like + {'obs_ok': True, 'comments': 'some comment'} + :return: dict + dictionary with stats + """ + if not obs_status_override: + obs_status_override = {} + + if telem is None: + telem = get_telemetry_by_observations( + star_obs, ignore_exceptions=True, as_table=False + ) + + if len(telem) != len(star_obs): + raise ValueError( + f"Length of telem ({len(telem)}) does not match length of star_obs ({len(star_obs)})." + ) + + # exclude star_obs that are in obs_status_override with status != 0 + excluded_obs = np.array( + [ + ( + (oi, ai) in obs_status_override + and obs_status_override[(oi, ai)]["status"] != 0 + ) + for oi, ai in star_obs[["mp_starcat_time", "agasc_id"]] + ] + ) + if np.any(excluded_obs): + logger.debug( + " Excluding observations flagged in obs-status table: " + f"{list(star_obs[excluded_obs]['obsid'])}" + ) + + included_obs = np.array( + [ + ( + (oi, ai) in obs_status_override + and obs_status_override[(oi, ai)]["status"] == 0 + ) + for oi, ai in star_obs[["mp_starcat_time", "agasc_id"]] + ] + ) + if np.any(included_obs): + logger.debug( + " Including observations marked OK in obs-status table: " + f"{list(star_obs[included_obs]['obsid'])}" + ) + + failures = [] + stats = [] + for i, obs in enumerate(star_obs): + oi, ai = obs["mp_starcat_time", "agasc_id"] + comment = "" + if (oi, ai) in obs_status_override: + status = obs_status_override[(oi, ai)] + logger.debug( + f" overriding status for (AGASC ID {ai}, starcat time {oi}): " + f"{status['status']}, {status['comments']}" + ) + comment = status["comments"] + + obs_telem = telem[i] + if "error_code" in obs_telem: + fail = obs_telem["error_code"] > 2 + obs_stat = get_obs_stats(obs, telem=[]) + obs_stat.update( + { + "obs_suspect": False, + "obs_fail": fail, + "excluded": excluded_obs[i], + "no_telem": True, + "comments": comment + if excluded_obs[i] + else f"Error: {obs_telem['msg']}.", + } + ) + stats.append(obs_stat) + if fail: + failures.append(obs_telem) + else: + obs_stat = get_obs_stats( + obs, telem={k: obs_telem[k] for k in obs_telem.colnames} + ) + obs_ok = included_obs[i] | (~excluded_obs[i] & obs_stat["obs_ok"]) + obs_stat.update( + { + "obs_ok": obs_ok, + "obs_suspect": not obs_ok and not excluded_obs[i], + "obs_fail": False, + "excluded": excluded_obs[i], + "no_telem": False, + "comments": comment, + } + ) + stats.append(obs_stat) + if obs_stat["obs_suspect"]: + failures.append( + dict( + MagStatsException( + msg="Suspect observation", + agasc_id=obs["agasc_id"], + obsid=obs["obsid"], + mp_starcat_time=obs["mp_starcat_time"], + ) + ) + ) + + return Table(stats), failures From da112c86d469be4db2cbc434652d9d3d4dfe529e Mon Sep 17 00:00:00 2001 From: Javier Gonzalez Date: Fri, 27 Feb 2026 13:28:23 -0500 Subject: [PATCH 03/14] WIP (a line that might not be needed) --- agasc/supplement/magnitudes/mag_estimate.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/agasc/supplement/magnitudes/mag_estimate.py b/agasc/supplement/magnitudes/mag_estimate.py index 23f41db..44c1fde 100644 --- a/agasc/supplement/magnitudes/mag_estimate.py +++ b/agasc/supplement/magnitudes/mag_estimate.py @@ -1622,6 +1622,8 @@ def get_multi_obs_stats(star_obs, telem=None, obs_status_override=None): if fail: failures.append(obs_telem) else: + # this conversion to Table might be needed to make sure different inputs work + obs_telem = Table(obs_telem) obs_stat = get_obs_stats( obs, telem={k: obs_telem[k] for k in obs_telem.colnames} ) From b13205bbf34d98708efb8f0a1bfe558037bc34f1 Mon Sep 17 00:00:00 2001 From: Javier Gonzalez Date: Tue, 17 Mar 2026 17:46:55 -0400 Subject: [PATCH 04/14] discard observations with no telemetry if they occured in the last two weeks --- agasc/supplement/magnitudes/mag_estimate.py | 28 ++++++++++++++++++--- 1 file changed, 25 insertions(+), 3 deletions(-) diff --git a/agasc/supplement/magnitudes/mag_estimate.py b/agasc/supplement/magnitudes/mag_estimate.py index 44c1fde..e00319b 100644 --- a/agasc/supplement/magnitudes/mag_estimate.py +++ b/agasc/supplement/magnitudes/mag_estimate.py @@ -17,6 +17,7 @@ from chandra_aca.transform import count_rate_to_mag, pixels_to_yagzag from cheta import fetch from cxotime import CxoTime +from cxotime import units as u from kadi import events from mica.archive import aca_l0 from mica.archive.aca_dark.dark_cal import get_dark_cal_image @@ -1178,7 +1179,9 @@ def calc_obs_stats(telem): } -def get_agasc_id_stats(agasc_id, obs_status_override=None, tstop=None): +def get_agasc_id_stats( + agasc_id, obs_status_override=None, tstop=None, no_telem_time=14 * u.day +): """ Get summary magnitude statistics for an AGASC ID. @@ -1189,6 +1192,8 @@ def get_agasc_id_stats(agasc_id, obs_status_override=None, tstop=None): {'obs_ok': True, 'comments': 'some comment'} :param tstop: cxotime-compatible timestamp Only entries in catalogs.STARS_OBS prior to this timestamp are considered. + :param no_telem_time: astropy Quantity + Time interval to consider for discarding recent observations with no telemetry. :return: dict dictionary with stats """ @@ -1204,11 +1209,28 @@ def get_agasc_id_stats(agasc_id, obs_status_override=None, tstop=None): ] star_obs.sort("mp_starcat_time") - n_obsids = len(star_obs) - all_telem = get_telemetry_by_observations( star_obs, ignore_exceptions=True, as_table=False ) + + # discard observations with no telemetry if they occured in the last two weeks. + discard = np.array( + [ + ("error_code" in tlm) + and (tlm["mp_starcat_time"] > CxoTime() - no_telem_time) + for tlm in all_telem + ] + ) + if np.any(discard): + obsids = [str(obsid) for obsid in star_obs["obsid"][discard]] + logger.debug( + f"Discarding recent observations with no telemetry (OBSIDs {' '.join(obsids)}." + ) + star_obs = star_obs[~discard] + all_telem = [tlm for sel, tlm in zip(~discard, all_telem, strict=True) if sel] + + n_obsids = len(star_obs) + stats, failures = get_multi_obs_stats( star_obs, obs_status_override=obs_status_override, telem=all_telem ) From c841850bcbd4a576f4e84a9b9ddd171a100ca872 Mon Sep 17 00:00:00 2001 From: Javier Gonzalez Date: Tue, 17 Mar 2026 17:47:49 -0400 Subject: [PATCH 05/14] remove catch for an exception that will not be raised --- agasc/supplement/magnitudes/update_mag_supplement.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/agasc/supplement/magnitudes/update_mag_supplement.py b/agasc/supplement/magnitudes/update_mag_supplement.py index 87aa92d..ed3cfc1 100755 --- a/agasc/supplement/magnitudes/update_mag_supplement.py +++ b/agasc/supplement/magnitudes/update_mag_supplement.py @@ -79,10 +79,6 @@ def get_agasc_id_stats( agasc_stats.append(agasc_stat) obs_stats.append(obs_stat) fails += obs_fail - except mag_estimate.MagStatsException as e: - msg = str(e) - logger.debug(msg) - fails.append(dict(e)) except Exception as e: # transform Exception to MagStatsException for standard book keeping msg = f"Unexpected Error: {e}" From ac4a0d4f5f93bc257d6825fef68f5c1e1abbbd01 Mon Sep 17 00:00:00 2001 From: Javier Gonzalez Date: Tue, 17 Mar 2026 18:09:37 -0400 Subject: [PATCH 06/14] added descriptions to sections in supplement report --- .../magnitudes/templates/run_report.html | 1 + .../magnitudes/update_mag_supplement.py | 20 +++++++++++++++++-- 2 files changed, 19 insertions(+), 2 deletions(-) diff --git a/agasc/supplement/magnitudes/templates/run_report.html b/agasc/supplement/magnitudes/templates/run_report.html index c2ca8a9..107a90a 100644 --- a/agasc/supplement/magnitudes/templates/run_report.html +++ b/agasc/supplement/magnitudes/templates/run_report.html @@ -56,6 +56,7 @@

{{ info.report_date }} Update Report

{%- for section in sections %}

{{ section.title }}

+

{{ section.description }}

diff --git a/agasc/supplement/magnitudes/update_mag_supplement.py b/agasc/supplement/magnitudes/update_mag_supplement.py index ed3cfc1..bc1a3da 100755 --- a/agasc/supplement/magnitudes/update_mag_supplement.py +++ b/agasc/supplement/magnitudes/update_mag_supplement.py @@ -729,10 +729,16 @@ def do( updated_stars["mag_aca_err"] != 0 ) sections = [ - {"id": "new_stars", "title": "New Stars", "stars": new_stars}, + { + "id": "new_stars", + "title": "New Stars", + "description": "These are stars that are being added to the supplement.", + "stars": new_stars, + }, { "id": "updated_stars", "title": "Updated Stars", + "description": "These are stars that are being updated in the supplement.", "stars": ( updated_stars["agasc_id"][updt_mag].tolist() if len(updated_stars[updt_mag]) @@ -742,6 +748,9 @@ def do( { "id": "not_updated_stars", "title": "Magnitude not Updated", + "description": ( + "These are stars whose magnitudes are not being updated in the supplement." + ), "stars": ( updated_stars["agasc_id"][~updt_mag].tolist() if len(updated_stars[~updt_mag]) @@ -750,7 +759,14 @@ def do( }, { "id": "other_stars", - "title": "Other (unexpectedly not updated)", + "title": "Stars in Limbo", + "description": ( + "These are stars that were in the list to process but are neither being" + " added or updated. This can happen if all observations for that star" + " fail or are skipped for some reason (e.g. a star with a single recent" + " observation that is suspect). This is resolved after the observations" + " are dispositioned or the failures are fixed." + ), "stars": list( agasc_stats["agasc_id"][ ~np.isin(agasc_stats["agasc_id"], new_stars) From cd5b779f5946487b49866c251f8997d636ea32dc Mon Sep 17 00:00:00 2001 From: Javier Gonzalez Date: Thu, 19 Mar 2026 10:01:43 -0400 Subject: [PATCH 07/14] make sure "excluded" and "no_telem" columns are present in obs_stats always, and then ignore them when writing the file --- agasc/supplement/magnitudes/mag_estimate.py | 3 +++ agasc/supplement/magnitudes/update_mag_supplement.py | 5 +++++ 2 files changed, 8 insertions(+) diff --git a/agasc/supplement/magnitudes/mag_estimate.py b/agasc/supplement/magnitudes/mag_estimate.py index e00319b..6cd439d 100644 --- a/agasc/supplement/magnitudes/mag_estimate.py +++ b/agasc/supplement/magnitudes/mag_estimate.py @@ -879,6 +879,8 @@ def get_obs_stats(obs, telem=None): "lf_variability_1000s": np.inf, "tempccd": np.nan, "dr_star": np.inf, + "no_telem": True, + "excluded": False, } ) @@ -973,6 +975,7 @@ def calc_obs_stats(telem): "n_mag_est_ok_3": n_mag_est_ok_3, "n_mag_est_ok_5": n_mag_est_ok_5, "dr_star": dr_star, + "no_telem": False, } if stats["n_mag_est_ok_3"] < 10: return stats diff --git a/agasc/supplement/magnitudes/update_mag_supplement.py b/agasc/supplement/magnitudes/update_mag_supplement.py index bc1a3da..e78c0a2 100755 --- a/agasc/supplement/magnitudes/update_mag_supplement.py +++ b/agasc/supplement/magnitudes/update_mag_supplement.py @@ -234,6 +234,11 @@ def update_mag_stats(obs_stats, agasc_stats, fails, outdir="."): if obs_stats is not None and len(obs_stats): filename = outdir / "mag_stats_obsid.fits" logger.debug(f"Updating {filename}") + + # these two were added late, and are not necessary in the file, so the file was not updated + columns = list(set(obs_stats.colnames) - {"excluded", "no_telem"}) + obs_stats = obs_stats[columns] + if filename.exists(): obs_stats = _update_table( table.Table.read(filename), obs_stats, keys=["agasc_id", "obsid"] From 7c91a60d90048323c03cfd6226162dd13f784850 Mon Sep 17 00:00:00 2001 From: Javier Gonzalez Date: Wed, 1 Apr 2026 15:23:48 -0400 Subject: [PATCH 08/14] in get_agasc_id_stats, handle the case where there is no valid stat out of get_multi_obs_stats --- agasc/supplement/magnitudes/mag_estimate.py | 23 +++++++++++++-------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/agasc/supplement/magnitudes/mag_estimate.py b/agasc/supplement/magnitudes/mag_estimate.py index 6cd439d..51c3c9c 100644 --- a/agasc/supplement/magnitudes/mag_estimate.py +++ b/agasc/supplement/magnitudes/mag_estimate.py @@ -1238,17 +1238,22 @@ def get_agasc_id_stats( star_obs, obs_status_override=obs_status_override, telem=all_telem ) - # combine magnitude estimates using a weighted mean - weighted_mean = get_weighted_mean(stats) - stats["w"] = weighted_mean["weights"] - stats["mean_corrected"] = weighted_mean["mean_corrected"] - stats["weighted_mean"] = weighted_mean["weighted_mean"] + # we do this in this method because a weighted mean doesn't make much sense for a list + # of observations. It only makes sense when considering all observations of a star. + if len(stats) > 0: + # combine magnitude estimates using a weighted mean + weighted_mean = get_weighted_mean(stats) + stats["w"] = weighted_mean["weights"] + stats["mean_corrected"] = weighted_mean["mean_corrected"] + stats["weighted_mean"] = weighted_mean["weighted_mean"] + else: + # or make sure column exists + stats["w"] = np.array([]) + stats["mean_corrected"] = np.array([]) + stats["weighted_mean"] = np.array([]) star = get_star(agasc_id, use_supplement=False) - # still need to check that this is the same as before - last_obs_time = CxoTime(stats["mp_starcat_time"][-1]).cxcsec - logger.debug(" identifying outlying observations...") for s, t in zip(stats, all_telem, strict=True): if s["no_telem"] or s["excluded"]: @@ -1333,7 +1338,7 @@ def get_agasc_id_stats( result.update( { "color": star["COLOR1"], - "last_obs_time": last_obs_time, + "last_obs_time": CxoTime(stats["mp_starcat_time"][-1]).cxcsec, "mag_aca": star["MAG_ACA"], "mag_aca_err": star["MAG_ACA_ERR"] / 100, "mag_obs_err": min_mag_obs_err, From 78459ea5013f6ce570f4f20faddc539a2d1b56a0 Mon Sep 17 00:00:00 2001 From: Javier Gonzalez Date: Wed, 1 Apr 2026 16:13:08 -0400 Subject: [PATCH 09/14] in get_telemetry_by_observations and get_telemetry_by_agasc_id, if as_table is True and no exception has been raised by the end, ignore all entries with errors, as this implies ignore_exceptions is true --- agasc/supplement/magnitudes/mag_estimate.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/agasc/supplement/magnitudes/mag_estimate.py b/agasc/supplement/magnitudes/mag_estimate.py index 51c3c9c..4c1208b 100644 --- a/agasc/supplement/magnitudes/mag_estimate.py +++ b/agasc/supplement/magnitudes/mag_estimate.py @@ -550,7 +550,13 @@ def get_telemetry_by_observations(observations, ignore_exceptions=False, as_tabl raise if telem: - return vstack([Table(tel) for tel in telem]) if as_table else telem + # we ignore entries with error_code here because this means ignore_exceptions is true + # (otherwise the exception would have been raised already) + return ( + vstack([Table(tel) for tel in telem if "error_code" not in tel]) + if as_table + else telem + ) return [] @@ -584,7 +590,8 @@ def add_obs_info(telem, obs_stats): o = telem["obsid"] == obsid telem["obs_ok"][o] = np.ones(np.count_nonzero(o), dtype=bool) * s["obs_ok"] if ( - np.any(telem["mag_est_ok"][o]) + len(telem[o]) > 0 + and np.any(telem["mag_est_ok"][o]) and s["f_mag_est_ok"] > 0 and np.isfinite(s["q75"]) and np.isfinite(s["q25"]) From 0ee4a5ac1c08a7842d32b0e8983dde940bdaddb0 Mon Sep 17 00:00:00 2001 From: Javier Gonzalez Date: Thu, 2 Apr 2026 10:40:43 -0400 Subject: [PATCH 10/14] Changed some comments to hopefully make things clearer in report. --- .../magnitudes/update_mag_supplement.py | 35 +++++++++++++++---- 1 file changed, 28 insertions(+), 7 deletions(-) diff --git a/agasc/supplement/magnitudes/update_mag_supplement.py b/agasc/supplement/magnitudes/update_mag_supplement.py index e78c0a2..c4a3a0d 100755 --- a/agasc/supplement/magnitudes/update_mag_supplement.py +++ b/agasc/supplement/magnitudes/update_mag_supplement.py @@ -260,8 +260,21 @@ def update_supplement(agasc_stats, filename, include_all=True, d_mag_threshold=0 """ Update the magnitude table of the AGASC supplement. - :param agasc_stats: - :param filename: + This function returns two list of stars: new stars and updated stars, which are stars that are + not in the supplement already, and stars thar are in the supplement before this update. Note + that "updated" stars are not necessarily updated in the supplement. If a star has a small change + in magnitude or magnitude uncertainty (less than d_mag_threshold), its magnitude is not updated + in the supplement, but it is still included in the "updated" stars list. + + This function compares the last_obs_time of the stars in the supplement and in agasc_stats. + If last_obs_time is the same, then the star is not updated in the supplement (and not included + in the "updated" stars list), even if there is a change in magnitude. This can happen if all + observations since last_updated are excluded or failed. + + :param agasc_stats: astropy.table.Table + The table with the new stats for each AGASC ID. It must include the columns in MAGS_DTYPE. + :param filename: str or pathlib.Path + The filename of the supplement to update. :param include_all: bool if True, all OK entries are included in supplement. if False, only OK entries marked 'selected_*' @@ -269,6 +282,8 @@ def update_supplement(agasc_stats, filename, include_all=True, d_mag_threshold=0 If the absolute difference between the new and the current mag_aca is less than this value, mag_aca is not updated. Note that last_obs_time is always updated. :return: + new_stars: list of AGASC IDs that are new in the supplement. + updated_stars: list of AGASC IDs that are already in the supplement. """ if agasc_stats is None or len(agasc_stats) == 0: return [], [] @@ -681,6 +696,10 @@ def do( except Exception as e: logger.warning(f"Failed to write {obs_status_file}: {e}") + # the following "updated_stars" is not the actual table of stars that were updated in the + # supplement, but the table of stars that would be updated if there were no threshold on d_mag. + # In other words: stars with tiny magnitude updates are not updated in the supplement, but they + # are included in "updated_stars". new_stars, updated_stars = update_supplement(agasc_stats, filename=filename) logger.info(f" {len(new_stars)} new stars, {len(updated_stars)} updated stars") @@ -766,11 +785,13 @@ def do( "id": "other_stars", "title": "Stars in Limbo", "description": ( - "These are stars that were in the list to process but are neither being" - " added or updated. This can happen if all observations for that star" - " fail or are skipped for some reason (e.g. a star with a single recent" - " observation that is suspect). This is resolved after the observations" - " are dispositioned or the failures are fixed." + "This section is here for informational purposes." + " These are stars that were in the list to process but are neither being" + " added nor updated. This can happen if all recent observations for that" + " star fail or are skipped for some reason (e.g. all recent observation are" + " so recent that telemetry is not available). This is resolved after the" + " observations are dispositioned or the failures are fixed." + " If there are errors, they should show up elsewhere in this report." ), "stars": list( agasc_stats["agasc_id"][ From 543515204d7c46fad3ce1b20988a1f001beb6584 Mon Sep 17 00:00:00 2001 From: Javier Gonzalez Date: Mon, 6 Apr 2026 15:44:55 -0400 Subject: [PATCH 11/14] Determine last mp_starcat_time with telemetry and cut processing right after that --- .../magnitudes/update_mag_supplement.py | 25 ++++++++++++++++++- 1 file changed, 24 insertions(+), 1 deletion(-) diff --git a/agasc/supplement/magnitudes/update_mag_supplement.py b/agasc/supplement/magnitudes/update_mag_supplement.py index c4a3a0d..b923b2c 100755 --- a/agasc/supplement/magnitudes/update_mag_supplement.py +++ b/agasc/supplement/magnitudes/update_mag_supplement.py @@ -566,6 +566,29 @@ def do( np.isin(star_obs_catalogs.STARS_OBS["agasc_id"], agasc_ids) ] + # find the latest observation with telemetry: take one star-obs per observation, + # and get telemetry for each just to find the last observation with data. + # Cut processing off right after that time. + recent_obs = stars_obs[stars_obs["mp_starcat_time"] > stop - 7 * u.day].copy() + recent_obs = recent_obs.group_by("mp_starcat_time") + recent_obs = recent_obs[recent_obs.groups.indices[:-1]] + recent_obs.sort(["mp_starcat_time"], reverse=True) + telem = mag_estimate.get_telemetry_by_observations( + recent_obs, ignore_exceptions=True, as_table=False + ) + processing_cutoff = stop + for obs, tel in zip(recent_obs, telem, strict=True): + if "error_code" in tel: + continue + logger.info( + f"Latest observation with telemetry: OBSID {obs['obsid']} at {obs['mp_starcat_time']}" + ) + processing_cutoff = CxoTime(obs["mp_starcat_time"]) + 1 * u.second + break + # and some AGASC Ids might be dropped because they are only observed in observations after the + # processing cutoff + agasc_ids = np.unique(agasc_ids) + # if supplement exists: # - drop bad stars # - get OBS status override @@ -662,7 +685,7 @@ def do( obs_stats, agasc_stats, fails = get_stats( agasc_ids, - tstop=stop, + tstop=processing_cutoff, obs_status_override=obs_status_override, no_progress=no_progress, ) From 80fe35525aa7fefde403203bd7922b63806d7652 Mon Sep 17 00:00:00 2001 From: Javier Gonzalez Date: Mon, 6 Apr 2026 16:02:12 -0400 Subject: [PATCH 12/14] fixup --- agasc/supplement/magnitudes/update_mag_supplement.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/agasc/supplement/magnitudes/update_mag_supplement.py b/agasc/supplement/magnitudes/update_mag_supplement.py index b923b2c..2bcdc2d 100755 --- a/agasc/supplement/magnitudes/update_mag_supplement.py +++ b/agasc/supplement/magnitudes/update_mag_supplement.py @@ -569,7 +569,9 @@ def do( # find the latest observation with telemetry: take one star-obs per observation, # and get telemetry for each just to find the last observation with data. # Cut processing off right after that time. - recent_obs = stars_obs[stars_obs["mp_starcat_time"] > stop - 7 * u.day].copy() + recent_obs = stars_obs[ + stars_obs["mp_starcat_time"] > CxoTime(stop) - 7 * u.day + ].copy() recent_obs = recent_obs.group_by("mp_starcat_time") recent_obs = recent_obs[recent_obs.groups.indices[:-1]] recent_obs.sort(["mp_starcat_time"], reverse=True) From 505020667fbe667d6651854a3cc57489128dda05 Mon Sep 17 00:00:00 2001 From: Javier Gonzalez Date: Wed, 15 Apr 2026 11:52:20 -0400 Subject: [PATCH 13/14] add a bit more logging --- agasc/supplement/magnitudes/update_mag_supplement.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/agasc/supplement/magnitudes/update_mag_supplement.py b/agasc/supplement/magnitudes/update_mag_supplement.py index 2bcdc2d..cb744ef 100755 --- a/agasc/supplement/magnitudes/update_mag_supplement.py +++ b/agasc/supplement/magnitudes/update_mag_supplement.py @@ -581,6 +581,9 @@ def do( processing_cutoff = stop for obs, tel in zip(recent_obs, telem, strict=True): if "error_code" in tel: + logger.info( + f"Skipping OBSID {obs['obsid']} at {obs['mp_starcat_time']} ({obs['msg']})" + ) continue logger.info( f"Latest observation with telemetry: OBSID {obs['obsid']} at {obs['mp_starcat_time']}" From 5ba8523fff4407bc09f836ad98c0f1cde1c43f3d Mon Sep 17 00:00:00 2001 From: Javier Gonzalez Date: Mon, 20 Apr 2026 09:55:00 -0400 Subject: [PATCH 14/14] fixup --- agasc/supplement/magnitudes/update_mag_supplement.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/agasc/supplement/magnitudes/update_mag_supplement.py b/agasc/supplement/magnitudes/update_mag_supplement.py index cb744ef..96b81ab 100755 --- a/agasc/supplement/magnitudes/update_mag_supplement.py +++ b/agasc/supplement/magnitudes/update_mag_supplement.py @@ -582,7 +582,8 @@ def do( for obs, tel in zip(recent_obs, telem, strict=True): if "error_code" in tel: logger.info( - f"Skipping OBSID {obs['obsid']} at {obs['mp_starcat_time']} ({obs['msg']})" + f"Skipping OBSID {obs['obsid']} at {obs['mp_starcat_time']}" + f" ({tel.get('msg', tel['error_code'])})" ) continue logger.info(