Chris773 commited on
Commit
3e1a13d
·
verified ·
1 Parent(s): 69d6c87

Update app.py

Browse files
Files changed (1) hide show
  1. app.py +198 -148
app.py CHANGED
@@ -19,17 +19,18 @@
19
  # TURN MODEL:
20
  # 3) Plausible final major turn (18:22→18:40) can adjust BOTH heading and speed
21
  #
22
- # SPEED MODEL (UPDATED):
23
  # 4) Southern trajectory speed is NOT constant:
24
  # - Solver chooses 2 speeds:
25
  # v1 from 18:40 → 19:40
26
  # v2 from 19:40 → 00:10
27
  # - Then 00:10 → 00:20 uses v2 (for display), but 00:20 is ignored in scoring.
28
  #
29
- # REQUESTED CHANGE (THIS UPDATE):
30
- # - Mirror the optimizer structure (continuous optimization with bounds) instead of beam search.
31
- # We now solve for: [south_heading, turn_speed, v1, v2]
32
- # using scipy.optimize.minimize (L-BFGS-B), with multi-start from a coarse grid.
 
33
  #
34
  # SCORING NOTE:
35
  # 5) 00:20 (24:20) is IGNORED for scoring:
@@ -51,6 +52,7 @@ import plotly.graph_objects as go
51
 
52
  import flight_plotter_v32 as fp
53
 
 
54
  from scipy.optimize import minimize
55
 
56
 
@@ -100,7 +102,6 @@ SATLOC_IDX_MEAS_BTO = 5
100
  RADAR_FIX_LAT = 6.578
101
  RADAR_FIX_LON = 96.341
102
  RADAR_FIX_TIME_UTC = datetime(2014, 3, 7, 18, 22, 0, tzinfo=timezone.utc)
103
-
104
  TURN_COMPLETE_TIME_UTC = datetime(2014, 3, 7, 18, 40, 0, tzinfo=timezone.utc)
105
 
106
  # Initial track shown in your radar annotation
@@ -118,7 +119,21 @@ GUIDED_VALIDATION_LABELS = [
118
  "24:20:00", # 00:20Z (final)
119
  ]
120
  GUIDED_FINAL_LABEL = "24:20:00" # 00:20Z FINAL RESTING POINT
121
- GUIDED_SCORE_EXCLUDE_LABELS = {GUIDED_FINAL_LABEL} # ignore 24:20 in scoring
 
 
 
 
 
 
 
 
 
 
 
 
 
 
122
 
123
  # Turn modeling parameters
124
  TURN_BANK_MIN_DEG = 15.0
@@ -126,13 +141,11 @@ TURN_BANK_MAX_DEG = 30.0
126
  TURN_BANK_DEFAULT_DEG = 25.0
127
  TURN_STEP_SEC = 30.0
128
 
129
- # Turn speed search parameters (UPDATED per request)
130
- TURN_SPEED_SEARCH_DELTA_KT = 40.0
131
- TURN_SPEED_STEP_KT = 5.0
132
  TURN_SPEED_MIN_ABS_KT = 470.0
133
  TURN_SPEED_MAX_ABS_KT = 500.0
134
 
135
- # Southern speed settings (2 segments + final display leg)
136
  SOUTH_SPEED_CENTER_GUESS_KT = 445.0
137
  SOUTH_SPEED_MIN_ABS_KT = 440.0
138
  SOUTH_SPEED_MAX_ABS_KT = 450.0
@@ -152,9 +165,16 @@ SPEED_CENTER_WEIGHT = 0.01
152
  GUIDED_HEADING_MIN = 150.0
153
  GUIDED_HEADING_MAX = 210.0
154
 
155
- # ===== Multi-start seed grid knobs (still cheap) =====
156
  COARSE_HEADING_STEP_DEG = 2.0
157
  COARSE_TURN_SPEED_STEP_KT = 5.0
 
 
 
 
 
 
 
158
 
159
  # Search area colors
160
  SEARCH_AREA_COLORS = {
@@ -671,7 +691,6 @@ def _map_marker_times(start_dt, legs):
671
  return times
672
 
673
 
674
- # ================= FAST PATH EVALUATOR (for optimizer) =================
675
  def _make_state_evaluator(start_lat: float, start_lon: float, legs, method: str, start_dt: datetime):
676
  """
677
  Precompute boundaries + leg start positions once, then evaluate lat/lon/spd/hdg at dt quickly.
@@ -1162,54 +1181,44 @@ def build_leaflet_map(path_pts, legs, method, start_dt):
1162
  """
1163
 
1164
 
1165
- # ================= SCORING HELPERS =================
1166
- def _build_pts_window_for_guided(start_dt: datetime, end_dt: datetime):
1167
- pts_all = _ping_points_after(start_dt)
1168
- pts_window = [(lbl, dt, bfo, bto) for (lbl, dt, bfo, bto) in pts_all if start_dt <= dt <= end_dt]
1169
- pts_window = [(lbl, dt, bfo, bto) for (lbl, dt, bfo, bto) in pts_window if lbl not in GUIDED_SCORE_EXCLUDE_LABELS]
1170
- return pts_window
1171
-
1172
-
1173
- def _regularization_speed(v1: float, v2: float) -> float:
1174
- # center + smooth penalties (kept light)
1175
- pen = 0.0
1176
- pen += float(SPEED_CENTER_WEIGHT) * ((float(v1) - float(SOUTH_SPEED_CENTER_GUESS_KT)) / 10.0) ** 2
1177
- pen += float(SPEED_CENTER_WEIGHT) * ((float(v2) - float(SOUTH_SPEED_CENTER_GUESS_KT)) / 10.0) ** 2
1178
- pen += float(SPEED_SMOOTH_WEIGHT) * ((float(v2) - float(v1)) / 10.0) ** 2
1179
- return float(pen)
1180
-
1181
-
1182
- def _compute_residuals_at_point(
1183
- method: str,
1184
- lat: float,
1185
- lon: float,
1186
- hdg_deg: float,
1187
- spd_knots: float,
1188
- dt_ping: datetime,
1189
- alt_ft: float,
1190
- vs_fpm: float,
1191
- descent_start_dt: datetime | None,
1192
- ):
1193
- alt_now_ft = altitude_ft_at_dt(dt_ping, RADAR_FIX_TIME_UTC, float(alt_ft), float(vs_fpm), descent_start_dt)
1194
- alt_m = alt_now_ft * 0.3048
1195
-
1196
- vs_mps = float(vs_fpm) * 0.00508
1197
- if descent_start_dt is not None and dt_ping < descent_start_dt:
1198
- vs_mps = 0.0
1199
-
1200
- spd_mps = float(spd_knots) * 1852.0 / 3600.0
1201
-
1202
- bfo_calc = float(fp.calculateBFO_at_time(spd_mps, float(hdg_deg), float(lat), float(lon), alt_m, vs_mps, dt_ping))
1203
- bto_calc = float(fp.calculateBTO_at_time(float(lat), float(lon), alt_m, dt_ping))
1204
-
1205
- bfo_meas, bto_meas = _lookup_measured_at_dt(dt_ping)
1206
-
1207
- r_bfo = None if bfo_meas is None else (bfo_calc - float(bfo_meas))
1208
- r_bto = None if bto_meas is None else (bto_calc - float(bto_meas))
1209
- return r_bfo, r_bto
1210
 
1211
 
1212
- # ================= Guided Solver (OPTIMIZER MIRROR) =================
1213
  def solve_guided_straight(
1214
  method, alt_ft, vs_fpm,
1215
  descent_start_choice, weight_bfo, weight_bto
@@ -1223,58 +1232,116 @@ def solve_guided_straight(
1223
  elif (descent_start_choice or "").startswith("After Ping Ring 7"):
1224
  descent_start_dt = dt_ping7_utc()
1225
 
1226
- pts_window = _build_pts_window_for_guided(start_dt, end_dt)
1227
- if len(pts_window) < 2:
1228
- return None, "Guided solve: not enough ping points to score.", None
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1229
 
1230
- # --- Objective: weighted normalized SSE across ping points (excluding 00:20) ---
1231
- def objective(x):
1232
- # x = [h_south, spd_turn, v1, v2]
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1233
  h_south = float(x[0])
1234
  spd_turn = float(x[1])
1235
  v1 = float(x[2])
1236
  v2 = float(x[3])
1237
 
1238
- # Build legs: turn (18:22->18:40) based on final heading & turn speed
1239
- turn_legs, _info = build_turn_legs(
1240
- speed_knots=spd_turn,
1241
- method=method,
1242
- heading_final_deg=h_south,
1243
- )
1244
- if not turn_legs:
1245
  return 1e12
1246
 
1247
- # Southern legs for scoring window only (18:40->00:10)
1248
- south_legs_score = [
1249
- (h_south, v1, float(SOUTH_SCHEDULE_SCORE_MIN[0])),
1250
- (h_south, v2, float(SOUTH_SCHEDULE_SCORE_MIN[1])),
1251
- ]
1252
-
1253
- # Composite legs up to 00:10 scoring horizon
1254
- legs_score = turn_legs + south_legs_score
1255
-
1256
- # Evaluator for states at ping times
1257
  state_at = _make_state_evaluator(RADAR_FIX_LAT, RADAR_FIX_LON, legs_score, method, start_dt)
1258
 
1259
- # Weighted normalized sum of squares
1260
  sse = 0.0
1261
  n_used = 0
1262
 
1263
- for (_lbl, dt_ping, _mbfo, _mbto) in pts_window:
1264
  st = state_at(dt_ping)
1265
  if st is None:
1266
  continue
1267
 
1268
- r_bfo, r_bto = _compute_residuals_at_point(
1269
- method=method,
1270
- lat=st["lat"],
1271
- lon=st["lon"],
1272
- hdg_deg=st["heading_deg"],
1273
- spd_knots=st["speed_kt"],
1274
- dt_ping=dt_ping,
1275
- alt_ft=alt_ft,
1276
- vs_fpm=vs_fpm,
1277
- descent_start_dt=descent_start_dt,
1278
  )
1279
 
1280
  if r_bfo is not None:
@@ -1287,38 +1354,45 @@ def solve_guided_straight(
1287
  if n_used == 0:
1288
  return 1e12
1289
 
1290
- # Light regularization (speed)
1291
  sse += _regularization_speed(v1, v2)
1292
-
1293
  return float(sse) / float(n_used)
1294
 
 
 
 
 
1295
  # Bounds (mirrors bounded optimizer style)
1296
  bounds = [
1297
- (GUIDED_HEADING_MIN, GUIDED_HEADING_MAX), # h_south
1298
- (TURN_SPEED_MIN_ABS_KT, TURN_SPEED_MAX_ABS_KT),# spd_turn
1299
- (SOUTH_SPEED_MIN_ABS_KT, SOUTH_SPEED_MAX_ABS_KT), # v1
1300
- (SOUTH_SPEED_MIN_ABS_KT, SOUTH_SPEED_MAX_ABS_KT), # v2
1301
  ]
1302
 
1303
- # Multi-start seeds (coarse grid only for initialization, then continuous optimize)
1304
  coarse_headings = np.arange(GUIDED_HEADING_MIN, GUIDED_HEADING_MAX + 1e-9, COARSE_HEADING_STEP_DEG)
1305
  coarse_turn_speeds = np.arange(TURN_SPEED_MIN_ABS_KT, TURN_SPEED_MAX_ABS_KT + 1e-9, COARSE_TURN_SPEED_STEP_KT)
1306
 
1307
- # Keep seed count reasonable
1308
- seeds = []
1309
  for h in coarse_headings:
1310
  for spd in coarse_turn_speeds:
1311
- seeds.append([float(h), float(spd), float(SOUTH_SPEED_CENTER_GUESS_KT), float(SOUTH_SPEED_CENTER_GUESS_KT)])
 
 
1312
 
1313
- best_res = None
 
1314
 
1315
- for x0 in seeds:
 
 
 
1316
  res = minimize(
1317
  objective,
1318
- x0=np.array(x0, dtype=float),
1319
  method="L-BFGS-B",
1320
  bounds=bounds,
1321
- options={"maxiter": 250, "ftol": 1e-9},
1322
  )
1323
  if not np.isfinite(res.fun):
1324
  continue
@@ -1333,31 +1407,13 @@ def solve_guided_straight(
1333
  v1_best = float(best_res.x[2])
1334
  v2_best = float(best_res.x[3])
1335
 
1336
- # Build final legs for display through 00:20 (00:20 excluded from scoring, but shown)
1337
- preturn_legs, _turn_info = build_turn_legs(
1338
- speed_knots=float(spd_turn_best),
1339
- method=method,
1340
- heading_final_deg=float(h_best),
1341
- )
1342
-
1343
- full_southern = [
1344
- (float(h_best), float(v1_best), float(SOUTH_SCHEDULE_FULL_MIN[0])),
1345
- (float(h_best), float(v2_best), float(SOUTH_SCHEDULE_FULL_MIN[1])),
1346
- (float(h_best), float(v2_best), float(SOUTH_SCHEDULE_FULL_MIN[2])),
1347
- ]
1348
-
1349
- full_legs = preturn_legs + full_southern
1350
  full_legs = clip_legs_to_end_dt(full_legs, start_dt, end_dt)
1351
 
1352
- # Compute final RMSEs for reporting (on scoring window)
1353
- # We evaluate on the same objective components, but report RMSE separately
1354
- # using raw residuals for BFO and BTO.
1355
- def final_metrics(h_south, spd_turn, v1, v2):
1356
- turn_legs, _ = build_turn_legs(spd_turn, method, h_south)
1357
- legs_score = turn_legs + [
1358
- (h_south, v1, float(SOUTH_SCHEDULE_SCORE_MIN[0])),
1359
- (h_south, v2, float(SOUTH_SCHEDULE_SCORE_MIN[1])),
1360
- ]
1361
  state_at = _make_state_evaluator(RADAR_FIX_LAT, RADAR_FIX_LON, legs_score, method, start_dt)
1362
 
1363
  sum_bfo2 = 0.0
@@ -1365,20 +1421,12 @@ def solve_guided_straight(
1365
  n_bfo = 0
1366
  n_bto = 0
1367
 
1368
- for (_lbl, dt_ping, _mbfo, _mbto) in pts_window:
1369
  st = state_at(dt_ping)
1370
  if st is None:
1371
  continue
1372
- r_bfo, r_bto = _compute_residuals_at_point(
1373
- method=method,
1374
- lat=st["lat"],
1375
- lon=st["lon"],
1376
- hdg_deg=st["heading_deg"],
1377
- spd_knots=st["speed_kt"],
1378
- dt_ping=dt_ping,
1379
- alt_ft=alt_ft,
1380
- vs_fpm=vs_fpm,
1381
- descent_start_dt=descent_start_dt,
1382
  )
1383
  if r_bfo is not None:
1384
  sum_bfo2 += float(r_bfo) ** 2
@@ -1391,12 +1439,13 @@ def solve_guided_straight(
1391
  rmse_bto = math.sqrt(sum_bto2 / max(1, n_bto)) if n_bto > 0 else None
1392
  return n_bfo, n_bto, rmse_bfo, rmse_bto
1393
 
1394
- n_bfo, n_bto, rmse_bfo, rmse_bto = final_metrics(h_best, spd_turn_best, v1_best, v2_best)
1395
 
1396
  report = "\n".join([
1397
  "**Guided mode: Straight + computed turn (18:22→18:40), clipped to 00:20**",
1398
- "- Optimizer: **scipy.optimize.minimize (L-BFGS-B)** with bounded variables (multi-start seeds)",
1399
- f"- Scoring excludes **00:20 (24:20)**",
 
1400
  f"- Best southern heading: **{h_best:.2f}°**",
1401
  f"- Best turn speed (18:22→18:40): **{spd_turn_best:.2f} kt**",
1402
  f"- Southern speed segments:",
@@ -1406,6 +1455,7 @@ def solve_guided_straight(
1406
  f"- Objective (normalized): **{float(best_res.fun):.6f}**",
1407
  f"- Optimizer status: **{best_res.message}**",
1408
  f"- Function evals: **{getattr(best_res, 'nfev', 'n/a')}** | Iters: **{getattr(best_res, 'nit', 'n/a')}**",
 
1409
  f"- Regularization: smooth={SPEED_SMOOTH_WEIGHT:.2f}, center={SPEED_CENTER_WEIGHT:.2f}",
1410
  f"- BFO points: **{n_bfo}** | RMSE_BFO: **{rmse_bfo:.2f} Hz**",
1411
  f"- BTO points: **{n_bto}** | RMSE_BTO: **{rmse_bto:.1f} µs**" if rmse_bto is not None else
@@ -1509,8 +1559,8 @@ with gr.Blocks() as demo:
1509
  f"- **Final resting point:** **00:20:00Z** (trajectory clipped exactly at this time)\n\n"
1510
  f"**Turn speed bounds:** **{TURN_SPEED_MIN_ABS_KT:.0f}–{TURN_SPEED_MAX_ABS_KT:.0f} kt**\n"
1511
  f"**Southern speed bounds:** **{SOUTH_SPEED_MIN_ABS_KT:.0f}–{SOUTH_SPEED_MAX_ABS_KT:.0f} kt** (center {SOUTH_SPEED_CENTER_GUESS_KT:.0f})\n"
1512
- f"**Scoring note:** **00:20 is excluded** from scoring and its measured values are removed from the validation table.\n"
1513
- f"**Solver:** bounded optimizer (L-BFGS-B) on [heading, turn speed, v1, v2] with coarse multi-start seeds."
1514
  )
1515
 
1516
  with gr.Row():
@@ -1530,7 +1580,7 @@ with gr.Blocks() as demo:
1530
  )
1531
 
1532
  gr.Markdown("### Guided mode (Straight only)")
1533
- gr.Markdown("- This build has **only** the straight guided solver (optimizer-based).")
1534
 
1535
  gr.Markdown("### Guided scoring weights")
1536
  with gr.Row():
 
19
  # TURN MODEL:
20
  # 3) Plausible final major turn (18:22→18:40) can adjust BOTH heading and speed
21
  #
22
+ # SPEED MODEL:
23
  # 4) Southern trajectory speed is NOT constant:
24
  # - Solver chooses 2 speeds:
25
  # v1 from 18:40 → 19:40
26
  # v2 from 19:40 → 00:10
27
  # - Then 00:10 → 00:20 uses v2 (for display), but 00:20 is ignored in scoring.
28
  #
29
+ # OPTIMIZER MIRROR + PERFORMANCE (THIS UPDATE):
30
+ # - Uses bounded optimizer (L-BFGS-B) on [heading, turn speed, v1, v2]
31
+ # - FAST seed pruning: pre-score coarse seeds using only 3 pings, keep top K
32
+ # - Strict scoring pings: only 18:40,19:40,20:40,21:40,22:40,00:10 (00:20 excluded)
33
+ # - Aggressive memoization cache for BFO/BTO calculations (big speed-up)
34
  #
35
  # SCORING NOTE:
36
  # 5) 00:20 (24:20) is IGNORED for scoring:
 
52
 
53
  import flight_plotter_v32 as fp
54
 
55
+ # Requires SciPy installed in your HF Space
56
  from scipy.optimize import minimize
57
 
58
 
 
102
  RADAR_FIX_LAT = 6.578
103
  RADAR_FIX_LON = 96.341
104
  RADAR_FIX_TIME_UTC = datetime(2014, 3, 7, 18, 22, 0, tzinfo=timezone.utc)
 
105
  TURN_COMPLETE_TIME_UTC = datetime(2014, 3, 7, 18, 40, 0, tzinfo=timezone.utc)
106
 
107
  # Initial track shown in your radar annotation
 
119
  "24:20:00", # 00:20Z (final)
120
  ]
121
  GUIDED_FINAL_LABEL = "24:20:00" # 00:20Z FINAL RESTING POINT
122
+
123
+ # Score ONLY these pings (and NEVER 00:20)
124
+ GUIDED_SCORE_LABELS = [
125
+ "18:40:00",
126
+ "19:40:00",
127
+ "20:40:00",
128
+ "21:40:00",
129
+ "22:40:00",
130
+ "24:10:00", # 00:10Z
131
+ ]
132
+ GUIDED_FAST_SEED_LABELS = [
133
+ "18:40:00",
134
+ "19:40:00",
135
+ "22:40:00",
136
+ ]
137
 
138
  # Turn modeling parameters
139
  TURN_BANK_MIN_DEG = 15.0
 
141
  TURN_BANK_DEFAULT_DEG = 25.0
142
  TURN_STEP_SEC = 30.0
143
 
144
+ # Turn speed bounds
 
 
145
  TURN_SPEED_MIN_ABS_KT = 470.0
146
  TURN_SPEED_MAX_ABS_KT = 500.0
147
 
148
+ # Southern speed bounds (2 segments)
149
  SOUTH_SPEED_CENTER_GUESS_KT = 445.0
150
  SOUTH_SPEED_MIN_ABS_KT = 440.0
151
  SOUTH_SPEED_MAX_ABS_KT = 450.0
 
165
  GUIDED_HEADING_MIN = 150.0
166
  GUIDED_HEADING_MAX = 210.0
167
 
168
+ # ===== Seed grid + pruning =====
169
  COARSE_HEADING_STEP_DEG = 2.0
170
  COARSE_TURN_SPEED_STEP_KT = 5.0
171
+ SEED_PRUNE_TOPK = 8 # BIG speed-up: only optimize from top K seeds
172
+
173
+ # Optimizer iterations
174
+ OPT_MAXITER = 160 # lower than before (faster), usually enough with good seeds
175
+
176
+ # Cache controls
177
+ CACHE_MAX_ITEMS = 50000 # safety cap
178
 
179
  # Search area colors
180
  SEARCH_AREA_COLORS = {
 
691
  return times
692
 
693
 
 
694
  def _make_state_evaluator(start_lat: float, start_lon: float, legs, method: str, start_dt: datetime):
695
  """
696
  Precompute boundaries + leg start positions once, then evaluate lat/lon/spd/hdg at dt quickly.
 
1181
  """
1182
 
1183
 
1184
+ # ================= Guided scoring point selection =================
1185
+ def _score_points_for_labels(labels: list[str]):
1186
+ """
1187
+ Return [(label, dt, meas_bfo, meas_bto)] using nearest measured sample for each label time.
1188
+ Ensures unique rows (if two labels snap to same satLoc index).
1189
+ """
1190
+ out = []
1191
+ used_idx = set()
1192
+ for lab in labels:
1193
+ dt = _dt_from_label(lab)
1194
+ if hasattr(fp, "nearest_measured_for_time"):
1195
+ idx = fp.nearest_measured_for_time(dt)[0]
1196
+ if idx in used_idx:
1197
+ continue
1198
+ used_idx.add(idx)
1199
+ row = fp.satLoc[idx]
1200
+ bfo = row[SATLOC_IDX_MEAS_BFO]
1201
+ bto = row[SATLOC_IDX_MEAS_BTO]
1202
+ if bfo is None:
1203
+ continue
1204
+ out.append((lab, dt, float(bfo), None if bto is None else float(bto)))
1205
+ else:
1206
+ # fallback: scan satLoc text label match
1207
+ best = None
1208
+ for r in fp.satLoc:
1209
+ label_r = fp._extract_time_string(r[0])
1210
+ if label_r.replace("Z", "") == lab.replace("Z", ""):
1211
+ best = r
1212
+ break
1213
+ if best is None or best[SATLOC_IDX_MEAS_BFO] is None:
1214
+ continue
1215
+ out.append((lab, dt, float(best[SATLOC_IDX_MEAS_BFO]),
1216
+ None if best[SATLOC_IDX_MEAS_BTO] is None else float(best[SATLOC_IDX_MEAS_BTO])))
1217
+ out.sort(key=lambda x: x[1])
1218
+ return out
 
 
 
 
 
 
 
 
 
 
1219
 
1220
 
1221
+ # ================= Guided Solver (OPTIMIZER MIRROR + FAST) =================
1222
  def solve_guided_straight(
1223
  method, alt_ft, vs_fpm,
1224
  descent_start_choice, weight_bfo, weight_bto
 
1232
  elif (descent_start_choice or "").startswith("After Ping Ring 7"):
1233
  descent_start_dt = dt_ping7_utc()
1234
 
1235
+ # Strict scoring pings (6 points)
1236
+ pts_score = _score_points_for_labels(GUIDED_SCORE_LABELS)
1237
+ if len(pts_score) < 2:
1238
+ return None, "Guided solve: not enough scoring pings (check satLoc).", None
1239
+
1240
+ pts_fast = _score_points_for_labels(GUIDED_FAST_SEED_LABELS)
1241
+ if len(pts_fast) < 2:
1242
+ pts_fast = pts_score[:3]
1243
+
1244
+ # --- Cache for expensive BFO/BTO computations ---
1245
+ cache = {}
1246
+
1247
+ def _cache_get(key):
1248
+ return cache.get(key, None)
1249
+
1250
+ def _cache_put(key, val):
1251
+ if len(cache) >= int(CACHE_MAX_ITEMS):
1252
+ # simple eviction: clear (safe + fast)
1253
+ cache.clear()
1254
+ cache[key] = val
1255
+
1256
+ def _compute_residuals_at_point_cached(lat, lon, hdg_deg, spd_knots, dt_ping):
1257
+ """
1258
+ Returns (r_bfo, r_bto) where residual = calc - meas.
1259
+ Uses memoization with rounding to speed repeated optimizer calls.
1260
+ """
1261
+ # altitude & VS depend on time
1262
+ alt_now_ft = altitude_ft_at_dt(dt_ping, RADAR_FIX_TIME_UTC, float(alt_ft), float(vs_fpm), descent_start_dt)
1263
+ alt_m = alt_now_ft * 0.3048
1264
 
1265
+ vs_mps = float(vs_fpm) * 0.00508
1266
+ if descent_start_dt is not None and dt_ping < descent_start_dt:
1267
+ vs_mps = 0.0
1268
+
1269
+ spd_mps = float(spd_knots) * 1852.0 / 3600.0
1270
+
1271
+ key = (
1272
+ round(float(lat), 4),
1273
+ round(float(lon), 4),
1274
+ round(float(hdg_deg), 1),
1275
+ round(float(spd_knots), 1),
1276
+ int(dt_ping.timestamp()),
1277
+ round(float(alt_m), 1),
1278
+ round(float(vs_mps), 3),
1279
+ )
1280
+ got = _cache_get(key)
1281
+ if got is not None:
1282
+ return got
1283
+
1284
+ bfo_calc = float(fp.calculateBFO_at_time(spd_mps, float(hdg_deg), float(lat), float(lon), alt_m, vs_mps, dt_ping))
1285
+ bto_calc = float(fp.calculateBTO_at_time(float(lat), float(lon), alt_m, dt_ping))
1286
+ bfo_meas, bto_meas = _lookup_measured_at_dt(dt_ping)
1287
+
1288
+ r_bfo = None if bfo_meas is None else (bfo_calc - float(bfo_meas))
1289
+ r_bto = None if bto_meas is None else (bto_calc - float(bto_meas))
1290
+
1291
+ _cache_put(key, (r_bfo, r_bto))
1292
+ return r_bfo, r_bto
1293
+
1294
+ def _regularization_speed(v1: float, v2: float) -> float:
1295
+ pen = 0.0
1296
+ pen += float(SPEED_CENTER_WEIGHT) * ((float(v1) - float(SOUTH_SPEED_CENTER_GUESS_KT)) / 10.0) ** 2
1297
+ pen += float(SPEED_CENTER_WEIGHT) * ((float(v2) - float(SOUTH_SPEED_CENTER_GUESS_KT)) / 10.0) ** 2
1298
+ pen += float(SPEED_SMOOTH_WEIGHT) * ((float(v2) - float(v1)) / 10.0) ** 2
1299
+ return float(pen)
1300
+
1301
+ def _build_legs_for_params(h_south, spd_turn, v1, v2, full=False):
1302
+ turn_legs, _ = build_turn_legs(
1303
+ speed_knots=float(spd_turn),
1304
+ method=method,
1305
+ heading_final_deg=float(h_south),
1306
+ )
1307
+ if not turn_legs:
1308
+ return None
1309
+ if full:
1310
+ south = [
1311
+ (float(h_south), float(v1), float(SOUTH_SCHEDULE_FULL_MIN[0])),
1312
+ (float(h_south), float(v2), float(SOUTH_SCHEDULE_FULL_MIN[1])),
1313
+ (float(h_south), float(v2), float(SOUTH_SCHEDULE_FULL_MIN[2])),
1314
+ ]
1315
+ else:
1316
+ south = [
1317
+ (float(h_south), float(v1), float(SOUTH_SCHEDULE_SCORE_MIN[0])),
1318
+ (float(h_south), float(v2), float(SOUTH_SCHEDULE_SCORE_MIN[1])),
1319
+ ]
1320
+ return turn_legs + south
1321
+
1322
+ # Objective on a chosen ping set
1323
+ def _objective_on_pts(x, pts):
1324
  h_south = float(x[0])
1325
  spd_turn = float(x[1])
1326
  v1 = float(x[2])
1327
  v2 = float(x[3])
1328
 
1329
+ legs_score = _build_legs_for_params(h_south, spd_turn, v1, v2, full=False)
1330
+ if legs_score is None:
 
 
 
 
 
1331
  return 1e12
1332
 
 
 
 
 
 
 
 
 
 
 
1333
  state_at = _make_state_evaluator(RADAR_FIX_LAT, RADAR_FIX_LON, legs_score, method, start_dt)
1334
 
 
1335
  sse = 0.0
1336
  n_used = 0
1337
 
1338
+ for (_lbl, dt_ping, _mbfo, _mbto) in pts:
1339
  st = state_at(dt_ping)
1340
  if st is None:
1341
  continue
1342
 
1343
+ r_bfo, r_bto = _compute_residuals_at_point_cached(
1344
+ st["lat"], st["lon"], st["heading_deg"], st["speed_kt"], dt_ping
 
 
 
 
 
 
 
 
1345
  )
1346
 
1347
  if r_bfo is not None:
 
1354
  if n_used == 0:
1355
  return 1e12
1356
 
 
1357
  sse += _regularization_speed(v1, v2)
 
1358
  return float(sse) / float(n_used)
1359
 
1360
+ # Full objective (strict scoring pings)
1361
+ def objective(x):
1362
+ return _objective_on_pts(x, pts_score)
1363
+
1364
  # Bounds (mirrors bounded optimizer style)
1365
  bounds = [
1366
+ (GUIDED_HEADING_MIN, GUIDED_HEADING_MAX), # heading
1367
+ (TURN_SPEED_MIN_ABS_KT, TURN_SPEED_MAX_ABS_KT), # turn speed
1368
+ (SOUTH_SPEED_MIN_ABS_KT, SOUTH_SPEED_MAX_ABS_KT), # v1
1369
+ (SOUTH_SPEED_MIN_ABS_KT, SOUTH_SPEED_MAX_ABS_KT), # v2
1370
  ]
1371
 
1372
+ # ===== Seed pruning (FAST) =====
1373
  coarse_headings = np.arange(GUIDED_HEADING_MIN, GUIDED_HEADING_MAX + 1e-9, COARSE_HEADING_STEP_DEG)
1374
  coarse_turn_speeds = np.arange(TURN_SPEED_MIN_ABS_KT, TURN_SPEED_MAX_ABS_KT + 1e-9, COARSE_TURN_SPEED_STEP_KT)
1375
 
1376
+ seed_scores = []
 
1377
  for h in coarse_headings:
1378
  for spd in coarse_turn_speeds:
1379
+ x0 = np.array([float(h), float(spd), float(SOUTH_SPEED_CENTER_GUESS_KT), float(SOUTH_SPEED_CENTER_GUESS_KT)], dtype=float)
1380
+ f0 = _objective_on_pts(x0, pts_fast)
1381
+ seed_scores.append((float(f0), float(h), float(spd)))
1382
 
1383
+ seed_scores.sort(key=lambda t: t[0])
1384
+ top = seed_scores[:max(1, int(SEED_PRUNE_TOPK))]
1385
 
1386
+ # ===== Optimize only top K seeds =====
1387
+ best_res = None
1388
+ for (_f0, h0, spd0) in top:
1389
+ x0 = np.array([float(h0), float(spd0), float(SOUTH_SPEED_CENTER_GUESS_KT), float(SOUTH_SPEED_CENTER_GUESS_KT)], dtype=float)
1390
  res = minimize(
1391
  objective,
1392
+ x0=x0,
1393
  method="L-BFGS-B",
1394
  bounds=bounds,
1395
+ options={"maxiter": int(OPT_MAXITER), "ftol": 1e-9},
1396
  )
1397
  if not np.isfinite(res.fun):
1398
  continue
 
1407
  v1_best = float(best_res.x[2])
1408
  v2_best = float(best_res.x[3])
1409
 
1410
+ # Build final legs for display through 00:20
1411
+ full_legs = _build_legs_for_params(h_best, spd_turn_best, v1_best, v2_best, full=True)
 
 
 
 
 
 
 
 
 
 
 
 
1412
  full_legs = clip_legs_to_end_dt(full_legs, start_dt, end_dt)
1413
 
1414
+ # Compute RMSEs for reporting (strict scoring pings)
1415
+ def final_metrics():
1416
+ legs_score = _build_legs_for_params(h_best, spd_turn_best, v1_best, v2_best, full=False)
 
 
 
 
 
 
1417
  state_at = _make_state_evaluator(RADAR_FIX_LAT, RADAR_FIX_LON, legs_score, method, start_dt)
1418
 
1419
  sum_bfo2 = 0.0
 
1421
  n_bfo = 0
1422
  n_bto = 0
1423
 
1424
+ for (_lbl, dt_ping, _mbfo, _mbto) in pts_score:
1425
  st = state_at(dt_ping)
1426
  if st is None:
1427
  continue
1428
+ r_bfo, r_bto = _compute_residuals_at_point_cached(
1429
+ st["lat"], st["lon"], st["heading_deg"], st["speed_kt"], dt_ping
 
 
 
 
 
 
 
 
1430
  )
1431
  if r_bfo is not None:
1432
  sum_bfo2 += float(r_bfo) ** 2
 
1439
  rmse_bto = math.sqrt(sum_bto2 / max(1, n_bto)) if n_bto > 0 else None
1440
  return n_bfo, n_bto, rmse_bfo, rmse_bto
1441
 
1442
+ n_bfo, n_bto, rmse_bfo, rmse_bto = final_metrics()
1443
 
1444
  report = "\n".join([
1445
  "**Guided mode: Straight + computed turn (18:22→18:40), clipped to 00:20**",
1446
+ "- Solver: **L-BFGS-B** on **[heading, turn_speed, v1, v2]**",
1447
+ f"- Seed pruning: coarse grid → **top {SEED_PRUNE_TOPK}** seeds (fast pre-score on {len(pts_fast)} pings)",
1448
+ f"- Strict scoring pings: **{', '.join(GUIDED_SCORE_LABELS)}**",
1449
  f"- Best southern heading: **{h_best:.2f}°**",
1450
  f"- Best turn speed (18:22→18:40): **{spd_turn_best:.2f} kt**",
1451
  f"- Southern speed segments:",
 
1455
  f"- Objective (normalized): **{float(best_res.fun):.6f}**",
1456
  f"- Optimizer status: **{best_res.message}**",
1457
  f"- Function evals: **{getattr(best_res, 'nfev', 'n/a')}** | Iters: **{getattr(best_res, 'nit', 'n/a')}**",
1458
+ f"- Cache size (end of solve): **{len(cache)}** items (cap {CACHE_MAX_ITEMS})",
1459
  f"- Regularization: smooth={SPEED_SMOOTH_WEIGHT:.2f}, center={SPEED_CENTER_WEIGHT:.2f}",
1460
  f"- BFO points: **{n_bfo}** | RMSE_BFO: **{rmse_bfo:.2f} Hz**",
1461
  f"- BTO points: **{n_bto}** | RMSE_BTO: **{rmse_bto:.1f} µs**" if rmse_bto is not None else
 
1559
  f"- **Final resting point:** **00:20:00Z** (trajectory clipped exactly at this time)\n\n"
1560
  f"**Turn speed bounds:** **{TURN_SPEED_MIN_ABS_KT:.0f}–{TURN_SPEED_MAX_ABS_KT:.0f} kt**\n"
1561
  f"**Southern speed bounds:** **{SOUTH_SPEED_MIN_ABS_KT:.0f}–{SOUTH_SPEED_MAX_ABS_KT:.0f} kt** (center {SOUTH_SPEED_CENTER_GUESS_KT:.0f})\n"
1562
+ f"**Scoring pings:** **{', '.join(GUIDED_SCORE_LABELS)}** (00:20 excluded)\n"
1563
+ f"**Performance:** seed_prune_topk={SEED_PRUNE_TOPK}, opt_maxiter={OPT_MAXITER}, cache_cap={CACHE_MAX_ITEMS}"
1564
  )
1565
 
1566
  with gr.Row():
 
1580
  )
1581
 
1582
  gr.Markdown("### Guided mode (Straight only)")
1583
+ gr.Markdown("- Optimizer-based (L-BFGS-B) with fast seed pruning + caching.")
1584
 
1585
  gr.Markdown("### Guided scoring weights")
1586
  with gr.Row():