Perform a system-level failure investigation

import numpy as np
rng = np.random.default_rng(61)
fs, T = 20, 80.0; n = int(fs * T); t = np.arange(n) / fs
# ---- synthetic flight: the vehicle flies along +x at 2 m/s; a featureless wall segment from t = 40 s degrades VO; the estimate then drifts and the controller tracks the (wrong) estimate
x_ref = 2.0 * t                                                         # reference route
drift_v = np.where(t > 40, 0.25 * (t - 40) / 10, 0.0).clip(max=0.5)     # VO velocity bias that develops after t = 40 s (m/s, capped)
est_err = np.cumsum(drift_v) / fs                                       # estimator position error grows from the biased velocity
x_est = x_ref + 0.05 * rng.standard_normal(n)                           # what the vehicle believes: the controller holds the estimate on the plan, so it looks fine
x_true = x_ref - est_err                                                # the controller tracks the estimate, so the true position moves away from the route
vo_features = np.where(t > 40, 18 + 4 * rng.standard_normal(n), 140 + 12 * rng.standard_normal(n))
inl_ratio = np.where(t > 40, 0.55 + 0.05 * rng.standard_normal(n), 0.92 + 0.02 * rng.standard_normal(n))
nis_vis = np.where(t > 40, 1 + 0.6 * np.clip(t - 40, 0, 15) ** 1.2 / 5, 1.0) * rng.chisquare(3, n) / 3
gnss_ref = x_true + 0.3 * rng.standard_normal(n)                        # an independent reference (for the investigation, e.g. a surveyed marker or an RTK log)
mode = np.where(t < 52.0, 'offboard', 'altitude'); contact = t[np.argmax(x_true > 112.0)] if (x_true > 112.0).any() else None
# ---- end of synthetic data: replace with your log's arrays

# 1) layers: what did the vehicle believe, want, and do?
track_err = x_est - x_ref; truth_err = x_true - x_ref; belief_err = x_est - x_true
print(f"max |estimate - plan| {np.abs(track_err).max():.2f} m (looks fine); max |true - plan| {np.abs(truth_err).max():.2f} m; max |estimate - truth| {np.abs(belief_err).max():.2f} m")
# 2) first deviation of each layer (cause precedes effect)
first = lambda cond: float(t[np.argmax(cond)]) if cond.any() else None
t_feat = first(vo_features < 60); t_inl = first(inl_ratio < 0.7); nis_smooth = np.convolve(nis_vis, np.ones(20) / 20, mode='same'); t_nis = first(nis_smooth > 2.0); t_belief = first(np.abs(belief_err) > 0.5); t_true = first(np.abs(truth_err) > 1.0)
print(f"first deviations: features {t_feat:.2f} s, inlier ratio {t_inl:.2f} s, vision NIS {t_nis:.2f} s, belief error > 0.5 m {t_belief:.2f} s, true position error > 1 m {t_true:.2f} s")
# 3) barriers: monitor, pilot
monitor_threshold_features = 3                                          # the (too low) threshold in the flown configuration
t_monitor = first(vo_features < monitor_threshold_features); t_pilot = 52.0
print(f"monitor (features < {monitor_threshold_features}) first trips at {t_monitor if t_monitor is None else round(t_monitor, 2)} s; pilot takeover at {t_pilot:.1f} s; contact at {contact:.1f} s")
det_delay = t_pilot - t_feat; excursion = 2.0 * det_delay + 4.0 / (2 * 3.0)
print(f"time from first anomaly to takeover {det_delay:.1f} s; excursion estimate {excursion:.1f} m")
# 4) what would the better barrier have done?
t_better = first(vo_features < 60); print(f"a features < 60 threshold with 3 consecutive frames would have flagged at about {t_better + 3 / fs:.2f} s, {t_pilot - (t_better + 3 / fs):.1f} s before the takeover")
assert abs(t_feat - 40.0) < 0.3 and t_feat < t_belief < t_true and np.abs(track_err).max() < 1.0 < np.abs(truth_err).max()
assert t_monitor is None or t_monitor > t_pilot - 1.0
assert contact is not None and contact > t_pilot and t_better + 3 / fs < t_pilot - 10