"""Step 5 (robustness): injection-recovery tests. (a) Noiseless, non-linear test of the covariance forecast (first version): For a random sample of analysed short-arc objects, a radial acceleration A1 is injected into the numerically integrated trajectory at the real observation epochs; a gravity-only orbit (6 parameters) is then fitted to these noiseless synthetic positions by Gauss-Newton iteration with the same weights. The post-fit chi^2 is the exact (non-linear) chi^2 penalty that a gravity-only fit would pay; the linear forecast predicts (A1 / sigma_A1)^2. (b) [revision] End-to-end test of the *measurement* with noise. For members of the tested samples T_1 and T_0.1, a radial acceleration (A/m = 1, 0.1 or 0) is injected into the trajectory, Gaussian noise equal to the assumed astrometric sigma is added at the real observation epochs, a gravity-only orbit is fitted iteratively to the synthetic astrometry (this stands in for the catalogue orbit), and A1 is then measured about that orbit exactly as in s03 (linearised 6 + 1 parameter fit with outlier rejection). Recorded: the recovered value, its pull, and whether the object would have been "cleared" (3-sigma upper bound below the injected threshold), i.e. the false-clearance rate of the search. (c) [revision] Recovery of a radial signal with the two-body partials on multi-year arcs (see main_twobody). Outputs: data/injection_test.csv, data/injection_measurement.csv, data/injection_twobody_partials.csv """ import os import sys import time import numpy as np import pandas as pd import orbfisher as of import s03_fisher as s3 from common import DATA, A1_PER_AMR, a2_yark_max N_SAMPLE = 60 SEED = 20261013 def model(state, t_ref, obs, ng): """Astrometric direction and 6-parameter design matrix (arcsec-normalised rows).""" tt = of.utc_to_tt(obs["mjd_utc"]) Y, nfev, ok = of.propagate(np.asarray(state, float), t_ref, tt, ng, True, 1e-10) ep, ev = of.earth_helio(tt) obspos = ep + obs["off"] r = Y[0:3].T; v = Y[3:6].T dist = np.linalg.norm(r - obspos, axis=1) dvec = r - v * (dist / of.C_AUD)[:, None] - obspos dist = np.linalg.norm(dvec, axis=1) dr = Y[6:60].T.reshape(-1, 6, 9)[:, 0:3, :] ra, dec, Ha, Hd = of._design_rows(dvec, dist, dr) return ra, dec, Ha[:, :6], Hd[:, :6] def inject(state0, t_ref, obs, a1, n_iter=8): """Returns (post-fit chi2 of gravity-only fit, pre-fit chi2, iterations used).""" ra_t, dec_t, _, _ = model(state0, t_ref, obs, (a1, 0.0, 0.0)) sig = obs["sig"] * of.ARCSEC x = np.array(state0, float) chi_prev = None chi0 = None for it in range(n_iter): ra, dec, Ha, Hd = model(x, t_ref, obs, (0.0, 0.0, 0.0)) da = ((ra_t - ra + np.pi) % (2 * np.pi) - np.pi) * np.cos(dec) / sig dd = (dec_t - dec) / sig chi = float(np.sum(da ** 2) + np.sum(dd ** 2)) if chi0 is None: chi0 = chi if chi_prev is not None and abs(chi_prev - chi) < 1e-4 * max(chi, 1e-12): break chi_prev = chi A = np.vstack([Ha / sig[:, None], Hd / sig[:, None]]) b = np.concatenate([da, dd]) sc = np.sqrt((A ** 2).sum(axis=0)) dx, *_ = np.linalg.lstsq(A / sc, b, rcond=1e-12) x = x + dx / sc return chi, chi0, it def synth_obs(state0, t_ref, obs, a1, rng, noise=1.0): """Synthetic astrometry of a trajectory with radial parameter a1 at the real epochs.""" ra_t, dec_t, _, _ = model(state0, t_ref, obs, (a1, 0.0, 0.0)) sig = obs["sig"] * of.ARCSEC o = dict(obs) o["ra"] = np.degrees(ra_t + noise * sig * rng.standard_normal(len(sig)) / np.cos(dec_t)) % 360.0 o["dec"] = np.degrees(dec_t + noise * sig * rng.standard_normal(len(sig))) return o def measure_injected(state0, t_ref, obs, a1, rng, noise=1.0): o = synth_obs(state0, t_ref, obs, a1, rng, noise) g = s3.gauss_newton(state0, t_ref, o, a1=0.0, solve_a1=False, n_iter=10) # "catalogue" orbit out = of.numeric_fisher(g["state"], t_ref, o, None, ng=(0.0, 0.0, 0.0)) m = of.fit_a1(out["Ha"], out["Hd"], out["res_ra"], out["res_dec"], o["sig"], groups=o.get("grp")) m["lin_A1"], m["lin_sig"], m["lin_rms7"] = m.get("m_A1", np.nan), m.get("m_sig", np.nan), m.get("m_rms7", np.nan) nl = s3.refine_nonlinear(g["state"], t_ref, o, m) # as in s03: non-linear refit of significant values m["nonlinear_refit"] = bool(nl.get("m_nl_ok", False)) if m["nonlinear_refit"]: m["m_A1"], m["m_sig"], m["m_rms7"] = nl["m_nl_A1"], nl["m_nl_sig"], nl["m_nl_rms7"] rms_g = float(np.sqrt(g["chi2"] / max(2 * g["n"] - 6, 1))) return m, rms_g def main_measurement(n1=80, n01=40, seed=SEED + 1): print("[5b] end-to-end injection-recovery of the measurement (with noise)") plan = s3.load_plan() fr = pd.read_csv(os.path.join(DATA, "fisher_results.csv"), dtype={"pdes": str}) df = plan.merge(fr, on="pdes") amr_sig = df.sig_A1_opt / A1_PER_AMR base = (df.model == "numeric") & df.ocat.isin(["short_inarc", "short_std"]) & np.isfinite(df.sig_A1) & (df.n_obs <= 400) \ & df.A1.isna() p1 = df[base & (amr_sig <= 1.0 / 6.0)].sample(n1, random_state=seed) p01 = df[base & (amr_sig <= 0.1 / 6.0) & ~df.pdes.isin(p1.pdes)] p01 = p01.sample(min(n01, len(p01)), random_state=seed) store = s3.ObsStore() hz = s3.load_hz() rng = np.random.default_rng(seed) rows = [] t0 = time.time() for pool, amrs in ((p1, (1.0, 0.0)), (p01, (1.0, 0.1, 0.0))): for row in pool.itertuples(index=False): obs = store.get(row.pdes, numbered=(row.kind == "an"), api=bool(row.post_freeze)) if obs is None: continue if row.ocat == "short_inarc": state0 = of.elements_to_state(row.a, row.e, row.i, row.om, row.w, row.ma); t_ref = float(row.epoch_mjd) else: state0, jd0 = s3.hz_state(hz, row.pdes, float(row.jd_hz)) if state0 is None: continue t_ref = jd0 - 2400000.5 for amr in amrs: try: m, rms_g = measure_injected(state0, t_ref, obs, amr * A1_PER_AMR, rng) except Exception as exc: # noqa continue if "m_A1" not in m: continue sg = m["m_sig"] * max(1.0, m["m_rms7"]) ul = float(of.upper_bound(m["m_A1"], sg)) / A1_PER_AMR ul_all = float(of.upper_bound(m["m_A1_all"], m["m_sig_all"] * max(1.0, m["m_rms7_all"]))) / A1_PER_AMR ul_lono = m.get("m_lono_ul", np.inf) / A1_PER_AMR # resolution rules as in s04_limits.add_sensitivity if m["m_A1"] / sg >= of.RESOLVE_Z: ul_all = np.inf if not m.get("m_lono_A1", np.inf) / m.get("m_lono_sig", 1.0) < of.RESOLVE_Z: ul_lono = np.inf rows.append(dict(pdes=row.pdes, sbin=row.sbin, arc_obs=row.arc_obs, n_obs=row.n_obs, amr_injected=amr, amr_sigma_forecast=row.sig_A1_opt / A1_PER_AMR, s_forecast=amr * A1_PER_AMR / row.sig_A1_opt, amr_recovered=m["m_A1"] / A1_PER_AMR, amr_recovered_sigma=sg / A1_PER_AMR, amr_ul3=ul, amr_ul3_all=ul_all, amr_lono_ul=ul_lono, amr_ul_best=min(ul, ul_all, ul_lono), amr_ul3_no_margin=(m["m_A1"] + 3 * sg) / A1_PER_AMR, nonlinear_refit=m["nonlinear_refit"], amr_recovered_linear=m["lin_A1"] / A1_PER_AMR, amr_ul3_linear_no_margin=(m["lin_A1"] + 3 * m["lin_sig"] * max(1.0, m["lin_rms7"])) / A1_PER_AMR, in_fq_sample_R0p7=bool(amr * A1_PER_AMR / row.sig_A1_opt >= 1.4 * np.sqrt(2 * row.n_obs)), in_fq_sample_R1=bool(amr * A1_PER_AMR / row.sig_A1_opt >= 2.0 * np.sqrt(2 * row.n_obs)), pull=(m["m_A1"] - amr * A1_PER_AMR) / sg, rms_gravity_only=rms_g, rms6=m["m_rms6"], rms7=m["m_rms7"], n_rejected=m["m_n_all"] - m["m_n"])) out = pd.DataFrame(rows) out.to_csv(os.path.join(DATA, "injection_measurement.csv"), index=False) for amr, g in out.groupby("amr_injected"): x = amr if amr > 0 else None msg = " injected A/m = %g: n = %d, pull median %.2f, 16-84%% %.2f..%.2f, max |pull| %.1f" % ( amr, len(g), g.pull.median(), g.pull.quantile(0.16), g.pull.quantile(0.84), g.pull.abs().max()) if x: z = g.amr_recovered / g.amr_recovered_sigma msg += "; falsely cleared by the primary fit: %d, by any of the three rules: %d (without the 0.75 margin: %d; linear step only, no margin: %d)" % ( int((g.amr_ul3 < x).sum()), int((g.amr_ul_best < x).sum()), int((g.amr_ul3_no_margin < x).sum()), int((g.amr_ul3_linear_no_margin < x).sum())) msg += "; non-linear refit used: %d; recovered/injected median %.3f (min %.3f, max %.3f); linear step only: min %.3f, max %.3f; detected at >= 5 sigma: %d" % ( int(g.nonlinear_refit.sum()), (g.amr_recovered / x).median(), (g.amr_recovered / x).min(), (g.amr_recovered / x).max(), (g.amr_recovered_linear / x).min(), (g.amr_recovered_linear / x).max(), int((z >= 5).sum())) msg += "; rms of the gravity-only fit: median %.2f, below 1.4 in %d, below 2.0 in %d" % ( g.rms_gravity_only.median(), int((g.rms_gravity_only < 1.4).sum()), int((g.rms_gravity_only < 2.0).sum())) for col, lim in (("in_fq_sample_R0p7", 1.4), ("in_fq_sample_R1", 2.0)): h = g[g[col]] msg += "; members of the fit-quality sample (%s): %d, smallest gravity-only rms %.2f (criterion expects >= %.1f)" % ( col[-5:].replace("_", " "), len(h), h.rms_gravity_only.min() if len(h) else np.nan, lim) else: msg += "; not cleared at X = 1 by the primary fit: %d, by any rule: %d" % ( int((g.amr_ul3 >= 1.0).sum()), int((g.amr_ul_best >= 1.0).sum())) print(msg) print(" (%.0f s)" % (time.time() - t0)) return out def main(n_sample=N_SAMPLE): print("[5] injection-recovery test") plan = s3.load_plan() fr = pd.read_csv(os.path.join(DATA, "fisher_results.csv"), dtype={"pdes": str}) df = plan.merge(fr, on="pdes") # short-arc objects analysed with the numeric model, where A/m = 1 m2/kg is forecast at 2-200 sigma s_pred = A1_PER_AMR / df.sig_A1 pool = df[(df.model == "numeric") & df.ocat.isin(["short_inarc", "short_std"]) & np.isfinite(df.sig_A1) & (s_pred > 2) & (s_pred < 200) & (df.n_obs <= 400)] samp = pool.sample(min(n_sample, len(pool)), random_state=SEED) store = s3.ObsStore() hz = s3.load_hz() rows = [] t0 = time.time() for row in samp.itertuples(index=False): obs = store.get(row.pdes, numbered=(row.kind == "an"), api=bool(row.post_freeze)) if obs is None: continue if row.ocat == "short_inarc": state0 = of.elements_to_state(row.a, row.e, row.i, row.om, row.w, row.ma); t_ref = float(row.epoch_mjd) else: state0, jd0 = s3.hz_state(hz, row.pdes, float(row.jd_hz)) if state0 is None: continue t_ref = jd0 - 2400000.5 for amr in (0.1, 1.0): a1 = amr * A1_PER_AMR try: chi, chi0, it = inject(state0, t_ref, obs, a1) except Exception as exc: # noqa continue rows.append(dict(pdes=row.pdes, sbin=row.sbin, arc_obs=row.arc_obs, n_obs=row.n_obs, amr_injected=amr, s_forecast=a1 / row.sig_A1, dchi2_forecast=(a1 / row.sig_A1) ** 2, dchi2_nonlinear=chi, chi2_before_refit=chi0, iterations=it)) out = pd.DataFrame(rows) out["ratio"] = np.sqrt(out.dchi2_nonlinear / out.dchi2_forecast) out.to_csv(os.path.join(DATA, "injection_test.csv"), index=False) for amr, g in out.groupby("amr_injected"): lr = np.log10(g.ratio) print(" A/m = %g: n = %d, sqrt(dchi2_nonlinear)/s_forecast median %.3f, 16-84%% %.3f-%.3f, min %.3f, max %.3f" % (amr, len(g), 10 ** lr.median(), 10 ** lr.quantile(0.16), 10 ** lr.quantile(0.84), g.ratio.min(), g.ratio.max())) print(" (%.0f s)" % (time.time() - t0)) return out def main_twobody(n=40): """(c) Efficiency of the two-body partials used for the measurement on multi-year arcs of objects > 100 m. For the 40 random multi-apparition objects of the two-body/numerical cross-check, a unit A1 signal is built with the *numerical* partials, a gravity-only fit with the numerical state partials is removed (as an orbit computer with a complete model would do), and A1 is then estimated from that post-fit residual with the *two-body* partials: (i) A1 alone, (ii) A1 with A2 solved (the primary multi-apparition measurement). The same is done for a unit A2 signal, to measure how much of a Yarkovsky drift leaks into A1.""" print("[5c] two-body partials: recovery of a radial signal on multi-year arcs") plan = s3.load_plan() store = s3.ObsStore() hz = s3.load_hz() sel = plan[plan.selected] ex = (sel.sbin.isin(["1-10 m", "10-100 m"]) | sel.A1.notna()) & ~sel.ocat.isin(["short_inarc", "short_std"]) chk = sel[(sel.ocat == "long_unnumbered") & ~ex].sort_values("u").head(n) rows = [] t0 = time.time() def lsq(A, y): sc = np.sqrt((A ** 2).sum(0)) x, *_ = np.linalg.lstsq(A / sc, y, rcond=1e-12) return x / sc for row in chk.itertuples(index=False): obs = store.get(row.pdes, api=bool(row.post_freeze)) if obs is None: continue st, jd0 = s3.hz_state(hz, row.pdes, float(np.median(obs["mjd_utc"])) + 2400000.5) if st is None: continue num = of.numeric_fisher(st, jd0 - 2400000.5, obs, None) kep = of.kepler_fisher(of.elements_to_state(row.a, row.e, row.i, row.om, row.w, row.ma), float(row.epoch_mjd), obs, None) sr = obs["sig"] * of.ARCSEC An = np.vstack([num["Ha"][:, :8] / sr[:, None], num["Hd"][:, :8] / sr[:, None]]) Ak = np.vstack([kep["Ha"][:, :8] / sr[:, None], kep["Hd"][:, :8] / sr[:, None]]) r1 = An[:, 6] - An[:, :6] @ lsq(An[:, :6], An[:, 6]) # post-fit residual of a unit A1 r2 = An[:, 7] - An[:, :6] @ lsq(An[:, :6], An[:, 7]) # post-fit residual of a unit A2 x1 = lsq(Ak[:, :7], r1) # two-body, A1 alone x1j = lsq(Ak, r1) # two-body, A1 with A2 x2 = lsq(Ak[:, :7], r2) x2j = lsq(Ak, r2) rows.append(dict(pdes=row.pdes, arc_obs=obs["mjd_utc"][-1] - obs["mjd_utc"][0], n_obs=len(sr), recovered_over_injected=float(x1[6]), recovered_over_injected_A2_solved=float(x1j[6]), A1_leak_per_unit_A2_alone=float(x2[6]), A1_leak_per_unit_A2_solved=float(x2j[6]), A2_recovered_over_injected=float(x2j[7]), sigma_twobody_over_numeric=of.marginal_sigma(Ak[:, :7].T @ Ak[:, :7], range(7), 6) / of.marginal_sigma(An[:, :7].T @ An[:, :7], range(7), 6), sigma_joint_twobody_over_numeric=of.marginal_sigma(Ak.T @ Ak, range(8), 6) / of.marginal_sigma(An.T @ An, range(8), 6))) out = pd.DataFrame(rows) out.to_csv(os.path.join(DATA, "injection_twobody_partials.csv"), index=False) for col, lab in (("recovered_over_injected", "A1 alone"), ("recovered_over_injected_A2_solved", "A1 with A2 solved"), ("A2_recovered_over_injected", "A2 (unit A2 injected)")): e = out[col] print(" %s, n = %d: recovered/injected median %.3f, 16-84%% %.3f-%.3f, min %.3f, max %.3f" % (lab, len(out), e.median(), e.quantile(0.16), e.quantile(0.84), e.min(), e.max())) a2max = float(a2_yark_max(100.0, 2000.0)) # largest Yarkovsky A2 of a 100 m body of 2000 kg/m3 for col, lab in (("A1_leak_per_unit_A2_alone", "A1 alone"), ("A1_leak_per_unit_A2_solved", "A1 with A2 solved")): e = out[col].abs() * a2max / A1_PER_AMR print(" leak of the maximum Yarkovsky A2 of a 100 m body (%.2e au/d2) into A/m [m2/kg], %s: median %.2e, max %.2e" % (a2max, lab, e.median(), e.max())) print(" (%.0f s)" % (time.time() - t0)) return out if __name__ == "__main__": n = [int(a.split("=")[1]) for a in sys.argv if a.startswith("--n=")] main(n[0] if n else N_SAMPLE) main_measurement() main_twobody()