Skip to content

Commit 0eba9e5

Browse files
Mark Stalzerclaude
andcommitted
Add per-LH innovation gate for reflection rejection (TE-PROC-039)
When light-outlier-threshold > 0, compute the pre-update RMS innovation for each lighthouse batch via a dry-run of map_light_data (no state change). If rms_innovation > threshold * light_residuals_all, skip that LH for the current sync cycle. Unlike the EWMA-history approach, this fires on the frame of the reflection itself rather than requiring multiple bad frames to accumulate. Disabled by default (threshold=0); a starting value of 5.0 is recommended. Adds 6 OutlierGate property tests to residual_cascade_props.c covering: disabled states (threshold=0, cold start), fire/no-fire on either side of the threshold, RMS correctness, and single-spike detection. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
1 parent 04303e1 commit 0eba9e5

4 files changed

Lines changed: 206 additions & 0 deletions

File tree

docs/specs/tracking-engine-specs.md

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -39,6 +39,7 @@ Prefix: **TE**
3939
- [x] **TE-PROC-036**: While the tracked object is stationary (IMU variance below `kalman-stationary-*` thresholds), the system shall apply a Zero Velocity Update pseudo-measurement forcing velocity and acceleration toward zero.
4040
- [x] **TE-PROC-037**: When tracking is lost (no valid measurements for the configured timeout period), the system shall reset the Kalman state covariance to reflect high uncertainty.
4141
- [x] **TE-PROC-038**: The system shall scale each lighthouse's observation noise covariance R by the ratio of that lighthouse's EWMA residual to the fleet mean residual, so that lighthouses with above-average residuals receive proportionally less Kalman weight.
42+
- [x] **TE-PROC-039**: Where `light-outlier-threshold` is > 0, the system shall compute the pre-update RMS innovation for each lighthouse batch and skip that lighthouse's Kalman update for the current sync cycle if the RMS exceeds `light-outlier-threshold × light_residuals_all`.
4243

4344
## Kalman Lighthouse Tracker
4445

src/survive_kalman_tracker.c

Lines changed: 25 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -46,6 +46,9 @@ STRUCT_CONFIG_SECTION(SurviveKalmanTracker)
4646
"Minimum observations to allow light data into the kalman filter", 16, t->light_required_obs)
4747

4848
STRUCT_CONFIG_ITEM("light-max-error", "Maximum error to integrate into lightcap", -1, t->lightcap_max_error)
49+
STRUCT_CONFIG_ITEM("light-outlier-threshold",
50+
"Per-LH RMS innovation gate: skip a LH batch if its pre-update RMS residual "
51+
"exceeds this multiple of light_residuals_all (0 = disabled)", 0, t->light_outlier_threshold)
4952
STRUCT_CONFIG_ITEM("kalman-light-variance", "Variance of raw light sensor readings", 1e-2, t->light_var)
5053
STRUCT_CONFIG_ITEM("obs-cov-scale", "Covariance matrix scaling for obs",
5154
1, t->obs_cov_scale)
@@ -447,6 +450,28 @@ void survive_kalman_tracker_integrate_saved_light(SurviveKalmanTracker *tracker,
447450
.tracker = tracker,
448451
};
449452

453+
// @spec TE-PROC-039
454+
// Per-LH innovation gate: compute pre-update RMS residual via a dry-run
455+
// of the measurement model (no state change). If it exceeds
456+
// light_outlier_threshold × light_residuals_all, the LH is producing
457+
// measurements too far from the current state estimate to be trusted
458+
// (typically a reflection or severe miscalibration) and is skipped for
459+
// this sync cycle.
460+
if (tracker->light_outlier_threshold > 0 && tracker->light_residuals_all > 0) {
461+
CN_CREATE_STACK_VEC(y_dry, cnt);
462+
bool dry_ok = map_light_data(&cbctx, &Z, &tracker->model.state, &y_dry, NULL);
463+
if (dry_ok) {
464+
const FLT *yv = cn_as_const_vector(&y_dry);
465+
FLT sq = 0;
466+
for (int i = 0; i < cnt; i++) sq += yv[i] * yv[i];
467+
FLT rms = FLT_SQRT(sq / cnt);
468+
if (rms > tracker->light_outlier_threshold * tracker->light_residuals_all) {
469+
tracker->stats.lightcap_model_dropped++;
470+
continue;
471+
}
472+
}
473+
}
474+
450475
SurviveObject *so = tracker->so;
451476

452477
// @spec TE-PROC-038

src/survive_kalman_tracker.h

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -140,6 +140,7 @@ typedef struct SurviveKalmanTracker {
140140
FLT Lightcap_R;
141141

142142
FLT lightcap_max_error;
143+
FLT light_outlier_threshold;
143144
int light_rampin_length;
144145
bool use_error_for_lh_pos;
145146

src/test_cases/residual_cascade_props.c

Lines changed: 179 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -964,3 +964,182 @@ TEST(AdaptiveR, PerLhEWMAConvergesIndependently) {
964964
}
965965
return 0;
966966
}
967+
968+
// ── 9. Per-LH Innovation Gate Properties (TE-PROC-039) ──────────
969+
//
970+
// The outlier gate in survive_kalman_tracker.c fires when:
971+
// rms_innovation > light_outlier_threshold * light_residuals_all
972+
//
973+
// These tests exercise the gate decision logic in isolation. They do not
974+
// call map_light_data (which requires the full Kalman machinery) — instead
975+
// they test the RMS computation and threshold comparison that wraps the
976+
// dry-run call.
977+
978+
// Reproduce the gate decision from survive_kalman_tracker.c:
979+
// gate fires when rms > threshold * mean_res (and threshold > 0, mean_res > 0)
980+
static bool outlier_gate_fires(double rms, double threshold, double mean_res) {
981+
if (threshold <= 0 || mean_res <= 0)
982+
return false;
983+
return rms > threshold * mean_res;
984+
}
985+
986+
// Reproduce the RMS computation:
987+
// rms = sqrt(sum(y[i]^2) / cnt)
988+
static double compute_rms(const double *y, int cnt) {
989+
double sq = 0;
990+
for (int i = 0; i < cnt; i++) sq += y[i] * y[i];
991+
return sqrt(sq / cnt);
992+
}
993+
994+
// Property 1: Gate is disabled when threshold <= 0 (default off)
995+
TEST(OutlierGate, DisabledWhenThresholdZero) {
996+
unsigned seed = (unsigned)time(NULL);
997+
srand(seed);
998+
999+
for (int trial = 0; trial < 10000; trial++) {
1000+
double rms = rand_range(0.0, 100.0);
1001+
double threshold = rand_range(-10.0, 0.0); // disabled range
1002+
double mean_res = rand_range(1e-6, 1.0);
1003+
1004+
if (outlier_gate_fires(rms, threshold, mean_res)) {
1005+
fprintf(stderr, "DisabledWhenThresholdZero FAILED (seed=%u, trial=%d): "
1006+
"rms=%.4f threshold=%.4f mean_res=%.4f\n", seed, trial, rms, threshold, mean_res);
1007+
return -1;
1008+
}
1009+
}
1010+
return 0;
1011+
}
1012+
1013+
// Property 2: Gate is disabled when mean_res <= 0 (cold start: no fleet baseline yet)
1014+
TEST(OutlierGate, DisabledOnColdStart) {
1015+
unsigned seed = (unsigned)time(NULL);
1016+
srand(seed);
1017+
1018+
for (int trial = 0; trial < 10000; trial++) {
1019+
double rms = rand_range(0.0, 100.0);
1020+
double threshold = rand_range(1.0, 10.0); // enabled
1021+
double mean_res = 0.0; // cold start
1022+
1023+
if (outlier_gate_fires(rms, threshold, mean_res)) {
1024+
fprintf(stderr, "DisabledOnColdStart FAILED (seed=%u, trial=%d): "
1025+
"rms=%.4f threshold=%.4f\n", seed, trial, rms, threshold);
1026+
return -1;
1027+
}
1028+
}
1029+
return 0;
1030+
}
1031+
1032+
// Property 3: Gate fires when rms > threshold * mean_res
1033+
TEST(OutlierGate, FiresAboveThreshold) {
1034+
unsigned seed = (unsigned)time(NULL);
1035+
srand(seed);
1036+
1037+
for (int trial = 0; trial < 10000; trial++) {
1038+
double threshold = rand_range(1.0, 10.0);
1039+
double mean_res = rand_range(1e-4, 0.1);
1040+
// rms strictly above the trigger point
1041+
double rms = threshold * mean_res + rand_range(1e-6, 0.5);
1042+
1043+
if (!outlier_gate_fires(rms, threshold, mean_res)) {
1044+
fprintf(stderr, "FiresAboveThreshold FAILED (seed=%u, trial=%d): "
1045+
"rms=%.6f threshold=%.4f mean_res=%.6f trigger=%.6f\n",
1046+
seed, trial, rms, threshold, mean_res, threshold * mean_res);
1047+
return -1;
1048+
}
1049+
}
1050+
return 0;
1051+
}
1052+
1053+
// Property 4: Gate does NOT fire when rms <= threshold * mean_res
1054+
TEST(OutlierGate, DoesNotFireBelowThreshold) {
1055+
unsigned seed = (unsigned)time(NULL);
1056+
srand(seed);
1057+
1058+
for (int trial = 0; trial < 10000; trial++) {
1059+
double threshold = rand_range(1.0, 10.0);
1060+
double mean_res = rand_range(1e-4, 0.1);
1061+
// rms strictly below the trigger point
1062+
double rms = rand_range(0.0, threshold * mean_res - 1e-8);
1063+
if (rms < 0) rms = 0;
1064+
1065+
if (outlier_gate_fires(rms, threshold, mean_res)) {
1066+
fprintf(stderr, "DoesNotFireBelowThreshold FAILED (seed=%u, trial=%d): "
1067+
"rms=%.6f threshold=%.4f mean_res=%.6f trigger=%.6f\n",
1068+
seed, trial, rms, threshold, mean_res, threshold * mean_res);
1069+
return -1;
1070+
}
1071+
}
1072+
return 0;
1073+
}
1074+
1075+
// Property 5: RMS is non-negative and zero only for all-zero innovations
1076+
TEST(OutlierGate, RmsNonNegativeAndZeroOnlyWhenAllZero) {
1077+
unsigned seed = (unsigned)time(NULL);
1078+
srand(seed);
1079+
1080+
// All-zero case
1081+
{
1082+
double y[8] = {0};
1083+
double rms = compute_rms(y, 8);
1084+
if (rms != 0.0) {
1085+
fprintf(stderr, "RmsNonNegativeAndZeroOnlyWhenAllZero FAILED: all-zero rms=%.10f\n", rms);
1086+
return -1;
1087+
}
1088+
}
1089+
1090+
for (int trial = 0; trial < 10000; trial++) {
1091+
int cnt = 1 + rand() % 32;
1092+
double y[32];
1093+
bool all_zero = true;
1094+
for (int i = 0; i < cnt; i++) {
1095+
y[i] = rand_range(-1.0, 1.0);
1096+
if (y[i] != 0.0) all_zero = false;
1097+
}
1098+
double rms = compute_rms(y, cnt);
1099+
if (rms < 0) {
1100+
fprintf(stderr, "RmsNonNegativeAndZeroOnlyWhenAllZero FAILED: rms=%.6f < 0\n", rms);
1101+
return -1;
1102+
}
1103+
if (!all_zero && rms == 0.0) {
1104+
fprintf(stderr, "RmsNonNegativeAndZeroOnlyWhenAllZero FAILED: non-zero y but rms=0\n");
1105+
return -1;
1106+
}
1107+
}
1108+
return 0;
1109+
}
1110+
1111+
// Property 6: Single large innovation (reflection scenario) triggers the gate.
1112+
// A reflection might produce one sensor innovation >> normal; compute_rms
1113+
// should be above a reasonable threshold even for a single outlier sensor.
1114+
TEST(OutlierGate, SingleLargeInnovationTriggers) {
1115+
unsigned seed = (unsigned)time(NULL);
1116+
srand(seed);
1117+
1118+
int triggered = 0;
1119+
1120+
for (int trial = 0; trial < 1000; trial++) {
1121+
int cnt = 8 + rand() % 24; // 8-32 sensors
1122+
double mean_res = rand_range(0.005, 0.02); // typical fleet residual
1123+
double threshold = 5.0;
1124+
1125+
// All sensors normal except one reflection spike
1126+
double y[32];
1127+
for (int i = 0; i < cnt; i++)
1128+
y[i] = rand_range(-mean_res, mean_res);
1129+
int spike_sensor = rand() % cnt;
1130+
// Reflection: 20-100× the normal noise
1131+
y[spike_sensor] = mean_res * rand_range(20.0, 100.0);
1132+
1133+
double rms = compute_rms(y, cnt);
1134+
if (outlier_gate_fires(rms, threshold, mean_res)) {
1135+
triggered++;
1136+
}
1137+
}
1138+
1139+
// A 20-100× spike on any single sensor in a batch of 8-32 should be detectable
1140+
if (triggered < 800) {
1141+
fprintf(stderr, "SingleLargeInnovationTriggers FAILED: only %d/1000 triggered\n", triggered);
1142+
return -1;
1143+
}
1144+
return 0;
1145+
}

0 commit comments

Comments
 (0)