Chris773 commited on
Commit
3401409
·
verified ·
1 Parent(s): 2ad6f71

Update app.py

Browse files
Files changed (1) hide show
  1. app.py +275 -407
app.py CHANGED
@@ -2,8 +2,8 @@
2
  # MH370 Flight Simulator — Leaflet map (Folium) + Ping Rings (KML)
3
  #
4
  # Guided scenarios (ACTIVE):
5
- # A) Straight heading (single southern heading)
6
- # B) Variable Heading and Speed (southern headings AND speeds vary by segment)
7
  #
8
  # Key behavior:
9
  # - Last radar spot FIXED: 18:22:00Z @ 6.578N, 96.341E
@@ -12,21 +12,15 @@
12
  # - Guided validation timestamps: 18:22, 18:40, 19:40, 20:40, 21:40, 22:40, 00:10, 00:20
13
  #
14
  # THIS VERSION:
15
- # - Guided solver is BTO-first "ring-fit"
16
  # - IMPORTANT: If a timestamp has no measured BTO (e.g. 18:40 in your satLoc),
17
  # the solver automatically does NOT enforce BTO there.
18
  # - Turn speed range: 470–520 kt
19
- # - Southern segment speed range: 430–510 kt (per segment)
20
  #
21
- # PERFORMANCE UPDATE (Top-K polish for Straight heading):
22
- # - The expensive “polish” step is NOT run for every candidate anymore.
23
- # - We search fast first, keep top-K candidates, then polish ONLY those.
24
- #
25
- # NEW (Variable Heading and Speed mode):
26
- # - After 18:40, BOTH heading and speed can vary per ping-to-ping segment.
27
- # - Uses a fast beam search over heading choices + per-segment speed fit (root/abs-min),
28
- # prioritizing BTO at measured ping rings.
29
- # - Heading changes are constrained (small window) to keep runtime reasonable.
30
 
31
  import os
32
  import re
@@ -74,13 +68,6 @@ RING_TIMES = {
74
  "Ping Ring 7": "24:10:00Z",
75
  }
76
 
77
- # Normalization (used for combining BFO + BTO into a single score)
78
- BFO_SCALE_HZ = 5.0
79
- BTO_SCALE_US = 50.0
80
-
81
- DEFAULT_WEIGHT_BFO = 1.0
82
- DEFAULT_WEIGHT_BTO = 2.0
83
-
84
  # satLoc indices
85
  SATLOC_IDX_MEAS_BFO = 4
86
  SATLOC_IDX_MEAS_BTO = 5
@@ -120,7 +107,7 @@ TURN_SPEED_MAX_KT = 520.0
120
  SOUTH_SPEED_MIN_KT = 430.0
121
  SOUTH_SPEED_MAX_KT = 510.0
122
 
123
- # Guided heading search (south-ish) — used for Straight heading mode
124
  GUIDED_HEADING_MIN = 175
125
  GUIDED_HEADING_MAX = 190
126
  COARSE_HEADING_STEP_DEG = 2.0
@@ -131,33 +118,25 @@ REFINE_HEADING_STEP_DEG = 0.5
131
  GUIDED_SOUTH_SEG_MIN = [60.0, 60.0, 60.0, 60.0, 90.0, 10.0]
132
  GUIDED_SOUTH_SEG_LABELS_END = ["19:40:00", "20:40:00", "21:40:00", "22:40:00", "24:10:00", "24:20:00"]
133
 
134
- # HARD BTO labels list (will be AUTO-FILTERED to those with measured BTO)
135
  HARD_BTO_LABELS = ["18:40:00", "19:40:00", "20:40:00", "21:40:00", "22:40:00", "24:10:00", "24:20:00"]
136
 
137
- # ===== Performance knobs: Top-K polish (Straight heading mode) =====
138
- TOPK_POLISH = 5 # polish only the best K candidates globally
139
- ENABLE_POLISH = True # keep polish, but only for top-K
140
- POLISH_PASSES = 3 # fewer passes = faster
141
- POLISH_WINDOW_KT = 6.0 # +/- window around current segment speed
142
- POLISH_STEP_KT = 2.0 # step for coordinate search
143
- POLISH_OBJ_EPS = 1e-3 # secondary weight for sum(abs)
144
-
145
- # ===== Variable Heading and Speed mode knobs (fast) =====
146
- VAR_HEADING_GLOBAL_MIN = 170.0
147
- VAR_HEADING_GLOBAL_MAX = 195.0
148
- VAR_TURN_HEADING_COARSE_MIN = 175.0
149
- VAR_TURN_HEADING_COARSE_MAX = 190.0
150
- VAR_TURN_HEADING_STEP = 2.0 # coarse choice for heading at 18:40
151
- VAR_TURN_SPEED_STEP = 10.0 # coarse scan for turn speed
152
- VAR_BEAM_WIDTH = 8 # keep small for speed
153
- VAR_HEADING_STEP = 1.5 # candidate headings per segment
154
- VAR_HEADING_WINDOW_DEG = 6.0 # +/- window around previous segment heading
155
- VAR_FIRST_WINDOW_DEG = 8.0 # first segment a bit wider
156
- VAR_MAX_CAND_PER_STATE = 9 # cap expansions
157
- VAR_ENABLE_LIGHT_REFINEMENT = True # optional second-pass with smaller step
158
- VAR_REFINE_STEP = 1.0
159
- VAR_REFINE_WINDOW = 3.0
160
- VAR_TOPK_FINAL = 10 # keep top N full trajectories before final scoring
161
 
162
  # Search area colors (optional)
163
  SEARCH_AREA_COLORS = {
@@ -700,7 +679,7 @@ def clip_legs_to_end_dt(legs, start_dt: datetime, end_dt: datetime):
700
  return out
701
 
702
 
703
- # ================= BFO ANALYSIS (RESTORED) =================
704
  def _ping_points_after(start_dt: datetime):
705
  pts = []
706
  for row in fp.satLoc:
@@ -1222,7 +1201,6 @@ def _fit_segment_speed(
1222
  e_lo = float(err(vmin))
1223
  e_hi = float(err(vmax))
1224
 
1225
- # If sign change: bisection root (fast + tight)
1226
  if (e_lo < 0.0 and e_hi > 0.0) or (e_lo > 0.0 and e_hi < 0.0):
1227
  a, b = float(vmin), float(vmax)
1228
  fa, fb = e_lo, e_hi
@@ -1240,8 +1218,7 @@ def _fit_segment_speed(
1240
  lat1, lon1 = _simulate_one_leg_end(method, lat0, lon0, hdg, v_best, minutes)
1241
  return float(v_best), float(err(v_best)), (lat1, lon1)
1242
 
1243
- # No sign change: minimize |err| on bounds
1244
- v_best, _ = _golden_min_abs(err, vmin, vmax, iters=22)
1245
  lat1, lon1 = _simulate_one_leg_end(method, lat0, lon0, hdg, v_best, minutes)
1246
  return float(v_best), float(err(v_best)), (lat1, lon1)
1247
 
@@ -1268,16 +1245,25 @@ def _compute_bto_errors_at_enforced(
1268
  return errs
1269
 
1270
 
1271
- def _build_legs_from_segment_lists(turn_legs, headings, speeds, minutes):
1272
- south_legs = []
 
 
 
 
 
 
 
 
 
1273
  for h, v, m in zip(headings, speeds, minutes):
1274
  if float(m) <= 0:
1275
  continue
1276
- south_legs.append((float(h), float(v), float(m)))
1277
- return turn_legs + south_legs
1278
 
1279
 
1280
- # ================= Straight heading mode: polish (unchanged logic) =================
1281
  def _polish_south_speeds_topk(
1282
  method: str,
1283
  start_dt: datetime,
@@ -1300,38 +1286,35 @@ def _polish_south_speeds_topk(
1300
 
1301
  def build_full(speeds_now: list[float]):
1302
  lat_cur, lon_cur = turn_end_lat, turn_end_lon
1303
- south_errs_local = []
1304
  for mins, end_lbl, vseg in zip(GUIDED_SOUTH_SEG_MIN, GUIDED_SOUTH_SEG_LABELS_END, speeds_now):
1305
  dt_end = _dt_from_label(end_lbl)
1306
  lat_cur, lon_cur = _simulate_one_leg_end(method, lat_cur, lon_cur, h_south, vseg, mins)
1307
  meas = _lookup_measured_bto(dt_end)
1308
  if meas is None:
1309
- south_errs_local.append(float("nan"))
1310
  else:
1311
  bto_here = _bto_calc(method, lat_cur, lon_cur, dt_end, alt_ft, vs_fpm, start_dt, descent_start_dt)
1312
- south_errs_local.append(float(bto_here - meas))
1313
 
1314
- full_legs = clip_legs_to_end_dt(
1315
- turn_legs + [(float(h_south), float(v), float(m)) for v, m in zip(speeds_now, GUIDED_SOUTH_SEG_MIN)],
1316
- start_dt,
1317
- end_dt,
1318
- )
1319
- return full_legs, south_errs_local
1320
 
1321
  def objective(full_legs):
1322
  bto_errs = _compute_bto_errors_at_enforced(method, start_dt, alt_ft, vs_fpm, descent_start_dt, full_legs, enforced_dts)
1323
  if not bto_errs:
1324
- return 1e9, []
1325
  abs_arr = np.abs(np.array(bto_errs, dtype=float))
1326
  max_abs = float(np.max(abs_arr))
1327
  sum_abs = float(np.sum(abs_arr))
1328
- return float(max_abs + POLISH_OBJ_EPS * sum_abs), bto_errs
1329
 
1330
- full_legs, south_errs = build_full(speeds)
1331
- best_obj, _ = objective(full_legs)
1332
  best_speeds = speeds[:]
1333
  best_full = full_legs
1334
- best_south_errs = south_errs
1335
 
1336
  for _p in range(int(passes)):
1337
  improved_any = False
@@ -1339,6 +1322,7 @@ def _polish_south_speeds_topk(
1339
  v0 = float(best_speeds[i])
1340
  v_lo = _clamp(v0 - float(window_kt), SOUTH_SPEED_MIN_KT, SOUTH_SPEED_MAX_KT)
1341
  v_hi = _clamp(v0 + float(window_kt), SOUTH_SPEED_MIN_KT, SOUTH_SPEED_MAX_KT)
 
1342
  cand = list(np.arange(v_lo, v_hi + 1e-9, float(step_kt)))
1343
  if not cand:
1344
  continue
@@ -1346,37 +1330,34 @@ def _polish_south_speeds_topk(
1346
  local_best_obj = best_obj
1347
  local_best_speeds = best_speeds[:]
1348
  local_best_full = best_full
1349
- local_best_south_errs = best_south_errs
1350
 
1351
  for v_try in cand:
1352
  sp = best_speeds[:]
1353
  sp[i] = float(v_try)
1354
- full_try, south_errs_try = build_full(sp)
1355
- obj_try, _ = objective(full_try)
1356
  if obj_try < local_best_obj - 1e-9:
1357
  local_best_obj = obj_try
1358
  local_best_speeds = sp
1359
  local_best_full = full_try
1360
- local_best_south_errs = south_errs_try
1361
 
1362
  if local_best_obj < best_obj - 1e-9:
1363
  best_obj = local_best_obj
1364
  best_speeds = local_best_speeds
1365
  best_full = local_best_full
1366
- best_south_errs = local_best_south_errs
1367
  improved_any = True
1368
 
1369
  if not improved_any:
1370
  break
1371
 
1372
- return best_speeds, best_full, best_south_errs
1373
 
1374
 
1375
- # ================= Guided solver: Straight heading =================
1376
- def solve_guided_straight_heading(
1377
- method, alt_ft, vs_fpm,
1378
- descent_start_choice, weight_bfo, weight_bto
1379
- ):
1380
  start_dt = _guided_start_dt()
1381
  end_dt = _dt_from_label(GUIDED_FINAL_LABEL)
1382
 
@@ -1394,9 +1375,10 @@ def solve_guided_straight_heading(
1394
  enforced_labels = _available_hard_bto_labels()
1395
  if not enforced_labels:
1396
  return None, "Guided solve: no measured BTO available at any enforced ring times.", None
 
1397
  enforced_dts = [_dt_from_label(lbl) for lbl in enforced_labels]
1398
 
1399
- def _evaluate(h_south: float, v_turn: float, do_polish: bool = False):
1400
  h_south = float(h_south)
1401
  v_turn = float(v_turn)
1402
 
@@ -1431,14 +1413,11 @@ def solve_guided_straight_heading(
1431
  south_errs.append(float(err_seg))
1432
  lat_cur, lon_cur = float(lat1), float(lon1)
1433
 
1434
- full_legs = clip_legs_to_end_dt(
1435
- turn_legs + [(float(h_south), float(v), float(m)) for v, m in zip(south_speeds, GUIDED_SOUTH_SEG_MIN)],
1436
- start_dt,
1437
- end_dt,
1438
- )
1439
 
1440
- if do_polish and ENABLE_POLISH:
1441
- polished_speeds, polished_full, polished_seg_errs = _polish_south_speeds_topk(
1442
  method=method,
1443
  start_dt=start_dt,
1444
  end_dt=end_dt,
@@ -1449,13 +1428,7 @@ def solve_guided_straight_heading(
1449
  turn_legs=turn_legs,
1450
  enforced_dts=enforced_dts,
1451
  init_speeds=south_speeds,
1452
- passes=POLISH_PASSES,
1453
- window_kt=POLISH_WINDOW_KT,
1454
- step_kt=POLISH_STEP_KT,
1455
  )
1456
- south_speeds = polished_speeds
1457
- full_legs = polished_full
1458
- south_errs = polished_seg_errs
1459
 
1460
  bto_errs = _compute_bto_errors_at_enforced(method, start_dt, alt_ft, vs_fpm, descent_start_dt, full_legs, enforced_dts)
1461
  if not bto_errs:
@@ -1467,10 +1440,8 @@ def solve_guided_straight_heading(
1467
  rmse = float(np.sqrt(np.mean(np.array(bto_errs, dtype=float) ** 2)))
1468
 
1469
  details = {
1470
- "mode": "Straight heading",
1471
  "h_south": float(h_south),
1472
  "v_turn": float(v_turn),
1473
- "south_headings": [float(h_south)] * len(GUIDED_SOUTH_SEG_MIN),
1474
  "south_speeds": [float(v) for v in south_speeds],
1475
  "south_errs_us": [float(e) for e in south_errs],
1476
  "turn_info": turn_info,
@@ -1478,18 +1449,15 @@ def solve_guided_straight_heading(
1478
  "max_abs_err_us": max_abs,
1479
  "sum_abs_err_us": sum_abs,
1480
  "rmse_bto_us": rmse,
1481
- "polished": bool(do_polish and ENABLE_POLISH),
1482
  }
1483
- key = (max_abs, sum_abs, rmse)
1484
- return key, full_legs, details
1485
 
1486
  def best_for_heading(h_south: float):
1487
- h_south = float(h_south)
1488
-
1489
  coarse_turns = list(np.arange(TURN_SPEED_MIN_KT, TURN_SPEED_MAX_KT + 1e-9, 5.0))
1490
  best_key = None
1491
  best_legs = None
1492
- best_details = None
1493
 
1494
  for v in coarse_turns:
1495
  out = _evaluate(h_south, v, do_polish=False)
@@ -1497,12 +1465,12 @@ def solve_guided_straight_heading(
1497
  continue
1498
  key, legs, det = out
1499
  if best_key is None or key < best_key:
1500
- best_key, best_legs, best_details = key, legs, det
1501
 
1502
  if best_key is None:
1503
  return None
1504
 
1505
- v0 = float(best_details["v_turn"])
1506
  v_lo = _clamp(v0 - 7.5, TURN_SPEED_MIN_KT, TURN_SPEED_MAX_KT)
1507
  v_hi = _clamp(v0 + 7.5, TURN_SPEED_MIN_KT, TURN_SPEED_MAX_KT)
1508
 
@@ -1513,32 +1481,33 @@ def solve_guided_straight_heading(
1513
  (max_abs, sum_abs, _rmse), _legs, _det = out2
1514
  return float(max_abs) + POLISH_OBJ_EPS * float(sum_abs)
1515
 
1516
- v_best, _ = _golden_min_abs(obj, v_lo, v_hi, iters=18)
1517
  out3 = _evaluate(h_south, v_best, do_polish=False)
1518
  if out3 is None:
1519
- return best_key, best_legs, best_details
1520
- key3, legs3, det3 = out3
1521
- return key3, legs3, det3
1522
 
 
1523
  coarse_headings = _frange(GUIDED_HEADING_MIN, GUIDED_HEADING_MAX, COARSE_HEADING_STEP_DEG)
1524
- global_best = None
1525
- global_legs = None
1526
- global_det = None
1527
-
1528
  candidates = []
 
 
 
 
1529
  for h in coarse_headings:
1530
  out = best_for_heading(h)
1531
  if out is None:
1532
  continue
1533
  key, legs, det = out
1534
  candidates.append((key, float(det["h_south"]), float(det["v_turn"])))
1535
- if global_best is None or key < global_best:
1536
- global_best, global_legs, global_det = key, legs, det
1537
 
1538
- if global_legs is None:
1539
- return None, "Guided solve: no solution found within speed bounds.", None
1540
 
1541
- h0 = float(global_det["h_south"])
 
1542
  refine_headings = _refine_range(h0, REFINE_HEADING_HALF_WINDOW_DEG, REFINE_HEADING_STEP_DEG, GUIDED_HEADING_MIN, GUIDED_HEADING_MAX)
1543
  for h in refine_headings:
1544
  out = best_for_heading(h)
@@ -1546,9 +1515,10 @@ def solve_guided_straight_heading(
1546
  continue
1547
  key, legs, det = out
1548
  candidates.append((key, float(det["h_south"]), float(det["v_turn"])))
1549
- if global_best is None or key < global_best:
1550
- global_best, global_legs, global_det = key, legs, det
1551
 
 
1552
  candidates_sorted = sorted(candidates, key=lambda x: x[0])
1553
  uniq = []
1554
  seen = set()
@@ -1561,20 +1531,15 @@ def solve_guided_straight_heading(
1561
  if len(uniq) >= max(1, int(TOPK_POLISH)):
1562
  break
1563
 
1564
- best_final_key = global_best
1565
- best_final_legs = global_legs
1566
- best_final_det = global_det
1567
-
1568
- if ENABLE_POLISH and uniq:
1569
  for _key, hh, vv in uniq:
1570
  outp = _evaluate(hh, vv, do_polish=True)
1571
  if outp is None:
1572
  continue
1573
  keyp, legsp, detp = outp
1574
  if best_final_key is None or keyp < best_final_key:
1575
- best_final_key = keyp
1576
- best_final_legs = legsp
1577
- best_final_det = detp
1578
 
1579
  south_speeds = best_final_det["south_speeds"]
1580
  south_errs = best_final_det["south_errs_us"]
@@ -1584,33 +1549,28 @@ def solve_guided_straight_heading(
1584
  "**Guided mode: Straight heading + computed turn (18:22→18:40), clipped to 00:20**",
1585
  "- Solver mode: **BTO-first ring-fit**",
1586
  f"- Best southern heading: **{best_final_det['h_south']:.2f}°**",
1587
- f"- Turn speed range: {TURN_SPEED_MIN_KT:.0f}–{TURN_SPEED_MAX_KT:.0f} kt | Best: **{best_final_det['v_turn']:.2f} kt**",
1588
- f"- Polish: **{'ON (Top-K only)' if ENABLE_POLISH else 'OFF'}** | Applied to winner: **{'YES' if polished_flag else 'NO'}**",
1589
  "",
1590
- "**Southern per-segment speeds (each constrained to 430–510 kt):**",
1591
- f"- 18:40→19:40: **{south_speeds[0]:.2f} kt** (segment end err {south_errs[0]:+.2f} µs)",
1592
- f"- 19:40→20:40: **{south_speeds[1]:.2f} kt** (segment end err {south_errs[1]:+.2f} µs)",
1593
- f"- 20:40→21:40: **{south_speeds[2]:.2f} kt** (segment end err {south_errs[2]:+.2f} µs)",
1594
- f"- 21:40→22:40: **{south_speeds[3]:.2f} kt** (segment end err {south_errs[3]:+.2f} µs)",
1595
- f"- 22:40→00:10: **{south_speeds[4]:.2f} kt** (segment end err {south_errs[4]:+.2f} µs)",
1596
- f"- 00:10→00:20: **{south_speeds[5]:.2f} kt** (segment end err {south_errs[5]:+.2f} µs)",
1597
  "",
1598
- f"- Enforced BTO timestamps (measured available): **{', '.join(best_final_det['enforced_labels'])}**",
1599
  f"- Max |BTO error| at enforced times: **{best_final_det['max_abs_err_us']:.2f} µs**",
1600
  f"- RMSE_BTO at enforced times: **{best_final_det['rmse_bto_us']:.2f} µs**",
1601
  "",
1602
  "*If 18:40 has no measured BTO in your satLoc, it is automatically not enforced.*",
1603
- f"*Performance: polish runs only for top-{TOPK_POLISH} candidates now.*",
1604
  ]
1605
-
1606
  return best_final_legs, "\n".join(report_lines), start_dt
1607
 
1608
 
1609
- # ================= Guided solver: Variable heading + speed (fast beam search) =================
1610
- def solve_guided_variable_heading_speed(
1611
- method, alt_ft, vs_fpm,
1612
- descent_start_choice, weight_bfo, weight_bto
1613
- ):
1614
  start_dt = _guided_start_dt()
1615
  end_dt = _dt_from_label(GUIDED_FINAL_LABEL)
1616
 
@@ -1625,239 +1585,174 @@ def solve_guided_variable_heading_speed(
1625
  return None, "Guided solve: no measured BTO available at any enforced ring times.", None
1626
  enforced_dts = [_dt_from_label(lbl) for lbl in enforced_labels]
1627
 
1628
- # -------- Beam search for south segments (after a given 18:40 heading) --------
1629
- def beam_south_from_turn_end(lat0, lon0, h0):
1630
- h0 = float(h0)
1631
-
1632
- def clamp_h(h):
1633
- return _clamp(float(h), VAR_HEADING_GLOBAL_MIN, VAR_HEADING_GLOBAL_MAX)
1634
-
1635
- # state: lat, lon, prev_h, headings[], speeds[], seg_end_errs[], key(max_abs,sum_abs)
1636
- beam = [{
1637
- "lat": float(lat0),
1638
- "lon": float(lon0),
1639
- "prev_h": float(h0),
1640
- "headings": [],
1641
- "speeds": [],
1642
- "seg_errs": [],
1643
- "max_abs": 0.0,
1644
- "sum_abs": 0.0,
1645
- }]
1646
-
1647
- for idx, (mins, end_lbl) in enumerate(zip(GUIDED_SOUTH_SEG_MIN, GUIDED_SOUTH_SEG_LABELS_END)):
1648
- dt_end = _dt_from_label(end_lbl)
1649
- next_states = []
1650
 
1651
- for st in beam:
1652
- prev_h = float(st["prev_h"])
1653
- window = VAR_FIRST_WINDOW_DEG if idx == 0 else VAR_HEADING_WINDOW_DEG
1654
-
1655
- # Candidate headings around prev_h
1656
- cand = _frange(clamp_h(prev_h - window), clamp_h(prev_h + window), VAR_HEADING_STEP)
1657
-
1658
- # If too many, subsample evenly (keeps speed bounded)
1659
- if len(cand) > VAR_MAX_CAND_PER_STATE:
1660
- take = VAR_MAX_CAND_PER_STATE
1661
- if take <= 1:
1662
- cand = [cand[len(cand)//2]]
1663
- else:
1664
- xs = np.linspace(0, len(cand) - 1, take)
1665
- cand = [cand[int(round(x))] for x in xs]
1666
-
1667
- for h_try in cand:
1668
- v_seg, err_seg, (lat1, lon1) = _fit_segment_speed(
1669
- method=method,
1670
- lat0=st["lat"],
1671
- lon0=st["lon"],
1672
- hdg=float(h_try),
1673
- minutes=float(mins),
1674
- dt_end=dt_end,
1675
- alt_ft=alt_ft,
1676
- vs_fpm=vs_fpm,
1677
- start_dt=start_dt,
1678
- descent_start_dt=descent_start_dt,
1679
- vmin=SOUTH_SPEED_MIN_KT,
1680
- vmax=SOUTH_SPEED_MAX_KT,
1681
- )
1682
-
1683
- # Use measured seg-end error if available; otherwise don't count it in partial score
1684
- meas = _lookup_measured_bto(dt_end)
1685
- if meas is None or (isinstance(err_seg, float) and np.isnan(err_seg)):
1686
- e_abs = 0.0
1687
- e_signed = float("nan")
1688
- else:
1689
- e_signed = float(err_seg)
1690
- e_abs = abs(float(err_seg))
1691
-
1692
- max_abs = max(float(st["max_abs"]), float(e_abs))
1693
- sum_abs = float(st["sum_abs"]) + float(e_abs)
1694
-
1695
- next_states.append({
1696
- "lat": float(lat1),
1697
- "lon": float(lon1),
1698
- "prev_h": float(h_try),
1699
- "headings": st["headings"] + [float(h_try)],
1700
- "speeds": st["speeds"] + [float(v_seg)],
1701
- "seg_errs": st["seg_errs"] + [float(e_signed) if not (isinstance(e_signed, float) and np.isnan(e_signed)) else float("nan")],
1702
- "max_abs": float(max_abs),
1703
- "sum_abs": float(sum_abs),
1704
- })
1705
-
1706
- if not next_states:
1707
- return []
1708
-
1709
- # Keep best by (max_abs, sum_abs) — segment-level heuristic
1710
- next_states.sort(key=lambda s: (s["max_abs"], s["sum_abs"]))
1711
- beam = next_states[:max(1, int(VAR_BEAM_WIDTH))]
1712
-
1713
- return beam
1714
-
1715
- # -------- Evaluate a candidate: choose 18:40 heading + turn speed, then beam search headings --------
1716
- def evaluate_candidate(h40, v_turn, do_refine=False):
1717
- h40 = float(h40)
1718
- v_turn = float(v_turn)
1719
 
1720
- turn_legs, turn_info = build_turn_legs(
1721
- speed_knots=float(v_turn),
1722
- method=method,
1723
- heading_final_deg=float(h40),
1724
- )
1725
- pts_turn, _, _ = simulate_path(RADAR_FIX_LAT, RADAR_FIX_LON, turn_legs, method)
1726
- lat_turn_end, lon_turn_end = float(pts_turn[-1][0]), float(pts_turn[-1][1])
1727
-
1728
- beams = beam_south_from_turn_end(lat_turn_end, lon_turn_end, h40)
1729
- if not beams:
1730
- return []
1731
-
1732
- finals = []
1733
- for st in beams:
1734
- full_legs = clip_legs_to_end_dt(
1735
- _build_legs_from_segment_lists(turn_legs, st["headings"], st["speeds"], GUIDED_SOUTH_SEG_MIN),
1736
- start_dt,
1737
- end_dt,
 
 
1738
  )
1739
- bto_errs = _compute_bto_errors_at_enforced(method, start_dt, alt_ft, vs_fpm, descent_start_dt, full_legs, enforced_dts)
1740
- if not bto_errs:
1741
- continue
1742
- abs_arr = np.abs(np.array(bto_errs, dtype=float))
1743
- max_abs = float(np.max(abs_arr))
1744
- sum_abs = float(np.sum(abs_arr))
1745
- rmse = float(np.sqrt(np.mean(np.array(bto_errs, dtype=float) ** 2)))
1746
-
1747
- finals.append({
1748
- "key": (max_abs, sum_abs, rmse),
1749
- "legs": full_legs,
1750
- "det": {
1751
- "mode": "Variable Heading and Speed",
1752
- "h40": float(h40),
1753
- "v_turn": float(v_turn),
1754
- "south_headings": [float(h) for h in st["headings"]],
1755
- "south_speeds": [float(v) for v in st["speeds"]],
1756
- "south_errs_us": [float(e) for e in st["seg_errs"]],
1757
- "turn_info": turn_info,
1758
- "enforced_labels": enforced_labels,
1759
- "max_abs_err_us": max_abs,
1760
- "sum_abs_err_us": sum_abs,
1761
- "rmse_bto_us": rmse,
1762
- }
1763
- })
1764
-
1765
- finals.sort(key=lambda x: x["key"])
1766
- finals = finals[:max(1, int(VAR_TOPK_FINAL))]
1767
-
1768
- # Optional refinement: re-run a tighter beam around the best headings (cheap-ish)
1769
- if do_refine and VAR_ENABLE_LIGHT_REFINEMENT and finals:
1770
- best = finals[0]["det"]
1771
- # Build tighter beams seeded by each segment heading (local windows)
1772
- # We re-run beam with smaller step/window, but still centered on previous heading (which follows best).
1773
- # To keep code simple, just temporarily override constants inside a nested closure.
1774
- old_step = VAR_HEADING_STEP
1775
- old_win = VAR_HEADING_WINDOW_DEG
1776
- old_first = VAR_FIRST_WINDOW_DEG
1777
- old_bw = VAR_BEAM_WIDTH
1778
- try:
1779
- # tighten
1780
- globals()["VAR_HEADING_STEP"] = float(VAR_REFINE_STEP)
1781
- globals()["VAR_HEADING_WINDOW_DEG"] = float(VAR_REFINE_WINDOW)
1782
- globals()["VAR_FIRST_WINDOW_DEG"] = float(VAR_REFINE_WINDOW + 1.0)
1783
- globals()["VAR_BEAM_WIDTH"] = max(6, int(VAR_BEAM_WIDTH)) # keep similar
1784
- # rerun from same turn end with the same h40/v_turn
1785
- finals2 = evaluate_candidate(h40, v_turn, do_refine=False)
1786
- if finals2:
1787
- finals = sorted(finals + finals2, key=lambda x: x["key"])[:max(1, int(VAR_TOPK_FINAL))]
1788
- finally:
1789
- globals()["VAR_HEADING_STEP"] = old_step
1790
- globals()["VAR_HEADING_WINDOW_DEG"] = old_win
1791
- globals()["VAR_FIRST_WINDOW_DEG"] = old_first
1792
- globals()["VAR_BEAM_WIDTH"] = old_bw
1793
-
1794
- return finals
1795
-
1796
- # -------- Coarse global search over (h40, v_turn) --------
1797
- h40_list = _frange(VAR_TURN_HEADING_COARSE_MIN, VAR_TURN_HEADING_COARSE_MAX, VAR_TURN_HEADING_STEP)
1798
- vturn_list = _frange(TURN_SPEED_MIN_KT, TURN_SPEED_MAX_KT, VAR_TURN_SPEED_STEP)
1799
-
1800
- best = None # dict with key/legs/det
1801
- pool = []
1802
-
1803
- for h40 in h40_list:
1804
- for vt in vturn_list:
1805
- finals = evaluate_candidate(h40, vt, do_refine=False)
1806
- if not finals:
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1807
  continue
1808
- # keep the best from this pair
1809
- pool.append(finals[0])
1810
- if best is None or finals[0]["key"] < best["key"]:
1811
- best = finals[0]
1812
-
1813
- if best is None:
1814
- return None, "Guided solve (Variable): no solution found within bounds.", None
1815
-
1816
- # Pick top few pairs and refine only those (fast)
1817
- pool.sort(key=lambda x: x["key"])
1818
- pool = pool[:5]
1819
- refined_best = best
1820
- for cand in pool:
1821
- det = cand["det"]
1822
- finals = evaluate_candidate(det["h40"], det["v_turn"], do_refine=True)
1823
- if finals and (finals[0]["key"] < refined_best["key"]):
1824
- refined_best = finals[0]
1825
-
1826
- det = refined_best["det"]
1827
- hs = det["south_headings"]
1828
- vs = det["south_speeds"]
1829
- es = det["south_errs_us"]
1830
-
1831
- def fmt_err(e):
1832
- if e is None:
1833
- return "n/a"
1834
- if isinstance(e, float) and np.isnan(e):
1835
- return "n/a"
1836
- return f"{float(e):+.2f} µs"
 
 
 
 
 
1837
 
1838
  report_lines = [
1839
- "**Guided mode: Variable Heading and Speed + computed turn (18:22→18:40), clipped to 00:20**",
1840
- "- Solver mode: **BTO-first ring-fit** (beam search headings + per-segment speed fit)",
1841
- f"- Heading @ 18:40 (post-turn): **{det['h40']:.2f}°**",
1842
- f"- Turn speed range: {TURN_SPEED_MIN_KT:.0f}–{TURN_SPEED_MAX_KT:.0f} kt | Best: **{det['v_turn']:.2f} kt**",
1843
- f"- Beam width: **{VAR_BEAM_WIDTH}** | Heading step: **{VAR_HEADING_STEP:.1f}°** | Window: **±{VAR_HEADING_WINDOW_DEG:.1f}°**",
1844
  "",
1845
- "**Per-segment heading + speed (end-of-segment BTO err where measured):**",
1846
- f"- 18:40→19:40: **{hs[0]:.1f}°**, **{vs[0]:.2f} kt** (err {fmt_err(es[0])})",
1847
- f"- 19:40→20:40: **{hs[1]:.1f}°**, **{vs[1]:.2f} kt** (err {fmt_err(es[1])})",
1848
- f"- 20:40→21:40: **{hs[2]:.1f}°**, **{vs[2]:.2f} kt** (err {fmt_err(es[2])})",
1849
- f"- 21:40→22:40: **{hs[3]:.1f}°**, **{vs[3]:.2f} kt** (err {fmt_err(es[3])})",
1850
- f"- 22:40→00:10: **{hs[4]:.1f}°**, **{vs[4]:.2f} kt** (err {fmt_err(es[4])})",
1851
- f"- 00:10→00:20: **{hs[5]:.1f}°**, **{vs[5]:.2f} kt** (err {fmt_err(es[5])})",
1852
  "",
1853
- f"- Enforced BTO timestamps (measured available): **{', '.join(det['enforced_labels'])}**",
1854
- f"- Max |BTO error| at enforced times: **{det['max_abs_err_us']:.2f} µs**",
1855
- f"- RMSE_BTO at enforced times: **{det['rmse_bto_us']:.2f} µs**",
1856
  "",
1857
  "*If 18:40 has no measured BTO in your satLoc, it is automatically not enforced.*",
1858
  ]
1859
-
1860
- return refined_best["legs"], "\n".join(report_lines), start_dt
1861
 
1862
 
1863
  # ================= RUNNERS =================
@@ -1899,23 +1794,11 @@ def run_sim_manual(start_lat, start_lon, legs_df, method, start_time_str, alt_ft
1899
  return gr.update(), "", html_map, stats_md, bfo_fig, bfo_md, alt_fig, gr.update(headers=headers, value=rows)
1900
 
1901
 
1902
- def run_guided(
1903
- guided_mode,
1904
- start_lat, start_lon, method, alt_ft, vs_fpm, descent_start_choice,
1905
- weight_bfo, weight_bto
1906
- ):
1907
- if (guided_mode or "").startswith("Variable"):
1908
- legs, report_md, start_dt = solve_guided_variable_heading_speed(
1909
- method,
1910
- alt_ft, vs_fpm,
1911
- descent_start_choice, weight_bfo, weight_bto
1912
- )
1913
  else:
1914
- legs, report_md, start_dt = solve_guided_straight_heading(
1915
- method,
1916
- alt_ft, vs_fpm,
1917
- descent_start_choice, weight_bfo, weight_bto
1918
- )
1919
 
1920
  if legs is None:
1921
  fig = go.Figure()
@@ -1961,9 +1844,15 @@ with gr.Blocks() as demo:
1961
  f"- **Initial track (fixed):** {RADAR_INITIAL_TRACK_DEG:.1f}°\n"
1962
  f"- **Final resting point:** **00:20:00Z** (trajectory clipped exactly at this time)\n\n"
1963
  f"**Turn speed bounds:** {TURN_SPEED_MIN_KT:.0f}–{TURN_SPEED_MAX_KT:.0f} kt\n"
1964
- f"**Southern speed bounds:** {SOUTH_SPEED_MIN_KT:.0f}–{SOUTH_SPEED_MAX_KT:.0f} kt (per segment)\n"
1965
- f"**Ring-fit goal:** match BTO at ping times where measured BTO exists (auto-detected from satLoc).\n"
1966
- f"**Straight heading perf:** polish runs only for **top-{TOPK_POLISH}** candidates."
 
 
 
 
 
 
1967
  )
1968
 
1969
  with gr.Row():
@@ -1982,23 +1871,6 @@ with gr.Blocks() as demo:
1982
  label="When does descent start? (VS applies after this time)",
1983
  )
1984
 
1985
- gr.Markdown("### Guided mode (Recommended)")
1986
- guided_mode = gr.Dropdown(
1987
- ["Straight heading", "Variable Heading and Speed"],
1988
- value="Straight heading",
1989
- label="Guided scenario",
1990
- )
1991
- gr.Markdown(
1992
- "- **Straight heading**: one southern heading after 18:40; speeds vary per segment; optional Top-K polish.\n"
1993
- "- **Variable Heading and Speed**: headings + speeds vary per segment; uses fast beam search to keep runtime reasonable.\n"
1994
- "- If 18:40 has no measured BTO in your satLoc, it is **automatically skipped** for enforcement."
1995
- )
1996
-
1997
- gr.Markdown("### Guided scoring weights (for reporting only)")
1998
- with gr.Row():
1999
- weight_bfo = gr.Slider(0.0, 5.0, value=DEFAULT_WEIGHT_BFO, step=0.1, label="Weight: BFO")
2000
- weight_bto = gr.Slider(0.0, 10.0, value=DEFAULT_WEIGHT_BTO, step=0.1, label="Weight: BTO")
2001
-
2002
  with gr.Row():
2003
  guided_btn = gr.Button("Generate: Guided scenario")
2004
  manual_btn = gr.Button("Simulate: Manual legs")
@@ -2045,11 +1917,7 @@ with gr.Blocks() as demo:
2045
 
2046
  guided_btn.click(
2047
  run_guided,
2048
- inputs=[
2049
- guided_mode,
2050
- start_lat, start_lon, method, alt_ft, vs_fpm, descent_start_choice,
2051
- weight_bfo, weight_bto
2052
- ],
2053
  outputs=[legs_df, status_md, out_map, stats, bfo_plot, bfo_score, alt_plot, validation_df],
2054
  )
2055
 
 
2
  # MH370 Flight Simulator — Leaflet map (Folium) + Ping Rings (KML)
3
  #
4
  # Guided scenarios (ACTIVE):
5
+ # A) Straight heading
6
+ # B) Variable heading & speed (beam search, BTO-first)
7
  #
8
  # Key behavior:
9
  # - Last radar spot FIXED: 18:22:00Z @ 6.578N, 96.341E
 
12
  # - Guided validation timestamps: 18:22, 18:40, 19:40, 20:40, 21:40, 22:40, 00:10, 00:20
13
  #
14
  # THIS VERSION:
15
+ # - Solver is BTO-first "ring-fit"
16
  # - IMPORTANT: If a timestamp has no measured BTO (e.g. 18:40 in your satLoc),
17
  # the solver automatically does NOT enforce BTO there.
18
  # - Turn speed range: 470–520 kt
19
+ # - Southern segment speed range: 430–510 kt
20
  #
21
+ # PERFORMANCE:
22
+ # - Straight heading mode keeps Top-K polish.
23
+ # - Variable heading & speed uses beam search and skips polish by default (much faster).
 
 
 
 
 
 
24
 
25
  import os
26
  import re
 
68
  "Ping Ring 7": "24:10:00Z",
69
  }
70
 
 
 
 
 
 
 
 
71
  # satLoc indices
72
  SATLOC_IDX_MEAS_BFO = 4
73
  SATLOC_IDX_MEAS_BTO = 5
 
107
  SOUTH_SPEED_MIN_KT = 430.0
108
  SOUTH_SPEED_MAX_KT = 510.0
109
 
110
+ # ===== Straight heading search =====
111
  GUIDED_HEADING_MIN = 175
112
  GUIDED_HEADING_MAX = 190
113
  COARSE_HEADING_STEP_DEG = 2.0
 
118
  GUIDED_SOUTH_SEG_MIN = [60.0, 60.0, 60.0, 60.0, 90.0, 10.0]
119
  GUIDED_SOUTH_SEG_LABELS_END = ["19:40:00", "20:40:00", "21:40:00", "22:40:00", "24:10:00", "24:20:00"]
120
 
121
+ # HARD BTO labels list (AUTO-FILTERED to those with measured BTO)
122
  HARD_BTO_LABELS = ["18:40:00", "19:40:00", "20:40:00", "21:40:00", "22:40:00", "24:10:00", "24:20:00"]
123
 
124
+ # ===== Performance knobs: Top-K polish (Straight mode only) =====
125
+ TOPK_POLISH = 5
126
+ ENABLE_POLISH_STRAIGHT = True
127
+ POLISH_PASSES = 3
128
+ POLISH_WINDOW_KT = 6.0
129
+ POLISH_STEP_KT = 2.0
130
+ POLISH_OBJ_EPS = 1e-3
131
+
132
+ # ===== Variable heading & speed solver knobs (FAST) =====
133
+ VAR_HEADING_MIN = 170.0
134
+ VAR_HEADING_MAX = 200.0
135
+ VAR_TURN_SPEED_STEP = 10.0 # coarse turn speed scan
136
+ VAR_FIRST_HEADING_STEP = 4.0 # coarse first south heading grid
137
+ VAR_HEADING_DELTA_SET = [-6, -3, 0, 3, 6] # candidate changes per segment (deg)
138
+ VAR_BEAM_WIDTH = 12 # keep top-N partial paths each segment
139
+ VAR_SMOOTH_PENALTY_US_PER_DEG = 1.5 # small penalty to avoid wiggles (still BTO-first)
 
 
 
 
 
 
 
 
140
 
141
  # Search area colors (optional)
142
  SEARCH_AREA_COLORS = {
 
679
  return out
680
 
681
 
682
+ # ================= BFO ANALYSIS =================
683
  def _ping_points_after(start_dt: datetime):
684
  pts = []
685
  for row in fp.satLoc:
 
1201
  e_lo = float(err(vmin))
1202
  e_hi = float(err(vmax))
1203
 
 
1204
  if (e_lo < 0.0 and e_hi > 0.0) or (e_lo > 0.0 and e_hi < 0.0):
1205
  a, b = float(vmin), float(vmax)
1206
  fa, fb = e_lo, e_hi
 
1218
  lat1, lon1 = _simulate_one_leg_end(method, lat0, lon0, hdg, v_best, minutes)
1219
  return float(v_best), float(err(v_best)), (lat1, lon1)
1220
 
1221
+ v_best, _ = _golden_min_abs(err, vmin, vmax, iters=18)
 
1222
  lat1, lon1 = _simulate_one_leg_end(method, lat0, lon0, hdg, v_best, minutes)
1223
  return float(v_best), float(err(v_best)), (lat1, lon1)
1224
 
 
1245
  return errs
1246
 
1247
 
1248
+ def _build_south_legs_piecewise_constant_heading(h_south: float, speeds: list[float], minutes: list[float]):
1249
+ legs = []
1250
+ for v, m in zip(speeds, minutes):
1251
+ if float(m) <= 0:
1252
+ continue
1253
+ legs.append((float(h_south), float(v), float(m)))
1254
+ return legs
1255
+
1256
+
1257
+ def _build_south_legs_variable(headings: list[float], speeds: list[float], minutes: list[float]):
1258
+ legs = []
1259
  for h, v, m in zip(headings, speeds, minutes):
1260
  if float(m) <= 0:
1261
  continue
1262
+ legs.append((float(h), float(v), float(m)))
1263
+ return legs
1264
 
1265
 
1266
+ # ================= Top-K polish (Straight mode only) =================
1267
  def _polish_south_speeds_topk(
1268
  method: str,
1269
  start_dt: datetime,
 
1286
 
1287
  def build_full(speeds_now: list[float]):
1288
  lat_cur, lon_cur = turn_end_lat, turn_end_lon
1289
+ seg_errs = []
1290
  for mins, end_lbl, vseg in zip(GUIDED_SOUTH_SEG_MIN, GUIDED_SOUTH_SEG_LABELS_END, speeds_now):
1291
  dt_end = _dt_from_label(end_lbl)
1292
  lat_cur, lon_cur = _simulate_one_leg_end(method, lat_cur, lon_cur, h_south, vseg, mins)
1293
  meas = _lookup_measured_bto(dt_end)
1294
  if meas is None:
1295
+ seg_errs.append(float("nan"))
1296
  else:
1297
  bto_here = _bto_calc(method, lat_cur, lon_cur, dt_end, alt_ft, vs_fpm, start_dt, descent_start_dt)
1298
+ seg_errs.append(float(bto_here - meas))
1299
 
1300
+ south_legs = _build_south_legs_piecewise_constant_heading(h_south, speeds_now, GUIDED_SOUTH_SEG_MIN)
1301
+ full_legs = clip_legs_to_end_dt(turn_legs + south_legs, start_dt, end_dt)
1302
+ return full_legs, seg_errs
 
 
 
1303
 
1304
  def objective(full_legs):
1305
  bto_errs = _compute_bto_errors_at_enforced(method, start_dt, alt_ft, vs_fpm, descent_start_dt, full_legs, enforced_dts)
1306
  if not bto_errs:
1307
+ return 1e9
1308
  abs_arr = np.abs(np.array(bto_errs, dtype=float))
1309
  max_abs = float(np.max(abs_arr))
1310
  sum_abs = float(np.sum(abs_arr))
1311
+ return float(max_abs + POLISH_OBJ_EPS * sum_abs)
1312
 
1313
+ full_legs, seg_errs = build_full(speeds)
1314
+ best_obj = objective(full_legs)
1315
  best_speeds = speeds[:]
1316
  best_full = full_legs
1317
+ best_seg_errs = seg_errs
1318
 
1319
  for _p in range(int(passes)):
1320
  improved_any = False
 
1322
  v0 = float(best_speeds[i])
1323
  v_lo = _clamp(v0 - float(window_kt), SOUTH_SPEED_MIN_KT, SOUTH_SPEED_MAX_KT)
1324
  v_hi = _clamp(v0 + float(window_kt), SOUTH_SPEED_MIN_KT, SOUTH_SPEED_MAX_KT)
1325
+
1326
  cand = list(np.arange(v_lo, v_hi + 1e-9, float(step_kt)))
1327
  if not cand:
1328
  continue
 
1330
  local_best_obj = best_obj
1331
  local_best_speeds = best_speeds[:]
1332
  local_best_full = best_full
1333
+ local_best_seg_errs = best_seg_errs
1334
 
1335
  for v_try in cand:
1336
  sp = best_speeds[:]
1337
  sp[i] = float(v_try)
1338
+ full_try, seg_try = build_full(sp)
1339
+ obj_try = objective(full_try)
1340
  if obj_try < local_best_obj - 1e-9:
1341
  local_best_obj = obj_try
1342
  local_best_speeds = sp
1343
  local_best_full = full_try
1344
+ local_best_seg_errs = seg_try
1345
 
1346
  if local_best_obj < best_obj - 1e-9:
1347
  best_obj = local_best_obj
1348
  best_speeds = local_best_speeds
1349
  best_full = local_best_full
1350
+ best_seg_errs = local_best_seg_errs
1351
  improved_any = True
1352
 
1353
  if not improved_any:
1354
  break
1355
 
1356
+ return best_speeds, best_full, best_seg_errs
1357
 
1358
 
1359
+ # ================= Guided Solver: Straight heading =================
1360
+ def solve_guided_straight(method, alt_ft, vs_fpm, descent_start_choice):
 
 
 
1361
  start_dt = _guided_start_dt()
1362
  end_dt = _dt_from_label(GUIDED_FINAL_LABEL)
1363
 
 
1375
  enforced_labels = _available_hard_bto_labels()
1376
  if not enforced_labels:
1377
  return None, "Guided solve: no measured BTO available at any enforced ring times.", None
1378
+
1379
  enforced_dts = [_dt_from_label(lbl) for lbl in enforced_labels]
1380
 
1381
+ def _evaluate(h_south: float, v_turn: float, do_polish: bool):
1382
  h_south = float(h_south)
1383
  v_turn = float(v_turn)
1384
 
 
1413
  south_errs.append(float(err_seg))
1414
  lat_cur, lon_cur = float(lat1), float(lon1)
1415
 
1416
+ south_legs = _build_south_legs_piecewise_constant_heading(float(h_south), south_speeds, GUIDED_SOUTH_SEG_MIN)
1417
+ full_legs = clip_legs_to_end_dt(turn_legs + south_legs, start_dt, end_dt)
 
 
 
1418
 
1419
+ if do_polish and ENABLE_POLISH_STRAIGHT:
1420
+ south_speeds, full_legs, south_errs = _polish_south_speeds_topk(
1421
  method=method,
1422
  start_dt=start_dt,
1423
  end_dt=end_dt,
 
1428
  turn_legs=turn_legs,
1429
  enforced_dts=enforced_dts,
1430
  init_speeds=south_speeds,
 
 
 
1431
  )
 
 
 
1432
 
1433
  bto_errs = _compute_bto_errors_at_enforced(method, start_dt, alt_ft, vs_fpm, descent_start_dt, full_legs, enforced_dts)
1434
  if not bto_errs:
 
1440
  rmse = float(np.sqrt(np.mean(np.array(bto_errs, dtype=float) ** 2)))
1441
 
1442
  details = {
 
1443
  "h_south": float(h_south),
1444
  "v_turn": float(v_turn),
 
1445
  "south_speeds": [float(v) for v in south_speeds],
1446
  "south_errs_us": [float(e) for e in south_errs],
1447
  "turn_info": turn_info,
 
1449
  "max_abs_err_us": max_abs,
1450
  "sum_abs_err_us": sum_abs,
1451
  "rmse_bto_us": rmse,
1452
+ "polished": bool(do_polish and ENABLE_POLISH_STRAIGHT),
1453
  }
1454
+ return (max_abs, sum_abs, rmse), full_legs, details
 
1455
 
1456
  def best_for_heading(h_south: float):
 
 
1457
  coarse_turns = list(np.arange(TURN_SPEED_MIN_KT, TURN_SPEED_MAX_KT + 1e-9, 5.0))
1458
  best_key = None
1459
  best_legs = None
1460
+ best_det = None
1461
 
1462
  for v in coarse_turns:
1463
  out = _evaluate(h_south, v, do_polish=False)
 
1465
  continue
1466
  key, legs, det = out
1467
  if best_key is None or key < best_key:
1468
+ best_key, best_legs, best_det = key, legs, det
1469
 
1470
  if best_key is None:
1471
  return None
1472
 
1473
+ v0 = float(best_det["v_turn"])
1474
  v_lo = _clamp(v0 - 7.5, TURN_SPEED_MIN_KT, TURN_SPEED_MAX_KT)
1475
  v_hi = _clamp(v0 + 7.5, TURN_SPEED_MIN_KT, TURN_SPEED_MAX_KT)
1476
 
 
1481
  (max_abs, sum_abs, _rmse), _legs, _det = out2
1482
  return float(max_abs) + POLISH_OBJ_EPS * float(sum_abs)
1483
 
1484
+ v_best, _ = _golden_min_abs(obj, v_lo, v_hi, iters=16)
1485
  out3 = _evaluate(h_south, v_best, do_polish=False)
1486
  if out3 is None:
1487
+ return best_key, best_legs, best_det
1488
+ return out3
 
1489
 
1490
+ # Stage 1
1491
  coarse_headings = _frange(GUIDED_HEADING_MIN, GUIDED_HEADING_MAX, COARSE_HEADING_STEP_DEG)
 
 
 
 
1492
  candidates = []
1493
+ best = None
1494
+ best_legs = None
1495
+ best_det = None
1496
+
1497
  for h in coarse_headings:
1498
  out = best_for_heading(h)
1499
  if out is None:
1500
  continue
1501
  key, legs, det = out
1502
  candidates.append((key, float(det["h_south"]), float(det["v_turn"])))
1503
+ if best is None or key < best:
1504
+ best, best_legs, best_det = key, legs, det
1505
 
1506
+ if best_legs is None:
1507
+ return None, "Guided solve: no solution found within bounds.", None
1508
 
1509
+ # Stage 2
1510
+ h0 = float(best_det["h_south"])
1511
  refine_headings = _refine_range(h0, REFINE_HEADING_HALF_WINDOW_DEG, REFINE_HEADING_STEP_DEG, GUIDED_HEADING_MIN, GUIDED_HEADING_MAX)
1512
  for h in refine_headings:
1513
  out = best_for_heading(h)
 
1515
  continue
1516
  key, legs, det = out
1517
  candidates.append((key, float(det["h_south"]), float(det["v_turn"])))
1518
+ if best is None or key < best:
1519
+ best, best_legs, best_det = key, legs, det
1520
 
1521
+ # Top-K polish
1522
  candidates_sorted = sorted(candidates, key=lambda x: x[0])
1523
  uniq = []
1524
  seen = set()
 
1531
  if len(uniq) >= max(1, int(TOPK_POLISH)):
1532
  break
1533
 
1534
+ best_final_key, best_final_legs, best_final_det = best, best_legs, best_det
1535
+ if ENABLE_POLISH_STRAIGHT and uniq:
 
 
 
1536
  for _key, hh, vv in uniq:
1537
  outp = _evaluate(hh, vv, do_polish=True)
1538
  if outp is None:
1539
  continue
1540
  keyp, legsp, detp = outp
1541
  if best_final_key is None or keyp < best_final_key:
1542
+ best_final_key, best_final_legs, best_final_det = keyp, legsp, detp
 
 
1543
 
1544
  south_speeds = best_final_det["south_speeds"]
1545
  south_errs = best_final_det["south_errs_us"]
 
1549
  "**Guided mode: Straight heading + computed turn (18:22→18:40), clipped to 00:20**",
1550
  "- Solver mode: **BTO-first ring-fit**",
1551
  f"- Best southern heading: **{best_final_det['h_south']:.2f}°**",
1552
+ f"- Turn speed: **{best_final_det['v_turn']:.2f} kt** (range {TURN_SPEED_MIN_KT:.0f}–{TURN_SPEED_MAX_KT:.0f})",
1553
+ f"- Polish: **{'ON (Top-K)' if ENABLE_POLISH_STRAIGHT else 'OFF'}** | Applied to winner: **{'YES' if polished_flag else 'NO'}**",
1554
  "",
1555
+ "**Southern per-segment speeds:**",
1556
+ f"- 18:40→19:40: **{south_speeds[0]:.2f} kt** (end err {south_errs[0]:+.2f} µs)",
1557
+ f"- 19:40→20:40: **{south_speeds[1]:.2f} kt** (end err {south_errs[1]:+.2f} µs)",
1558
+ f"- 20:40→21:40: **{south_speeds[2]:.2f} kt** (end err {south_errs[2]:+.2f} µs)",
1559
+ f"- 21:40→22:40: **{south_speeds[3]:.2f} kt** (end err {south_errs[3]:+.2f} µs)",
1560
+ f"- 22:40→00:10: **{south_speeds[4]:.2f} kt** (end err {south_errs[4]:+.2f} µs)",
1561
+ f"- 00:10→00:20: **{south_speeds[5]:.2f} kt** (end err {south_errs[5]:+.2f} µs)",
1562
  "",
1563
+ f"- Enforced BTO timestamps: **{', '.join(best_final_det['enforced_labels'])}**",
1564
  f"- Max |BTO error| at enforced times: **{best_final_det['max_abs_err_us']:.2f} µs**",
1565
  f"- RMSE_BTO at enforced times: **{best_final_det['rmse_bto_us']:.2f} µs**",
1566
  "",
1567
  "*If 18:40 has no measured BTO in your satLoc, it is automatically not enforced.*",
 
1568
  ]
 
1569
  return best_final_legs, "\n".join(report_lines), start_dt
1570
 
1571
 
1572
+ # ================= Guided Solver: Variable heading & speed (FAST beam search) =================
1573
+ def solve_guided_variable(method, alt_ft, vs_fpm, descent_start_choice):
 
 
 
1574
  start_dt = _guided_start_dt()
1575
  end_dt = _dt_from_label(GUIDED_FINAL_LABEL)
1576
 
 
1585
  return None, "Guided solve: no measured BTO available at any enforced ring times.", None
1586
  enforced_dts = [_dt_from_label(lbl) for lbl in enforced_labels]
1587
 
1588
+ # turn speed candidates (coarse)
1589
+ turn_speeds = list(np.arange(TURN_SPEED_MIN_KT, TURN_SPEED_MAX_KT + 1e-9, VAR_TURN_SPEED_STEP))
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1590
 
1591
+ # initial (first south heading) candidates (coarse)
1592
+ first_headings = _frange(VAR_HEADING_MIN, VAR_HEADING_MAX, VAR_FIRST_HEADING_STEP)
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1593
 
1594
+ # Beam state:
1595
+ # {
1596
+ # lat, lon, headings[], speeds[], v_turn, max_abs, sum_abs, smooth_sum
1597
+ # }
1598
+ best_global = None
1599
+ best_global_legs = None
1600
+ best_global_info = None
1601
+
1602
+ def score_from_errors(errs: list[float]):
1603
+ if not errs:
1604
+ return None
1605
+ abs_arr = np.abs(np.array(errs, dtype=float))
1606
+ return float(np.max(abs_arr)), float(np.sum(abs_arr)), float(np.sqrt(np.mean(abs_arr * abs_arr)))
1607
+
1608
+ for v_turn in turn_speeds:
1609
+ for h1 in first_headings:
1610
+ turn_legs, turn_info = build_turn_legs(
1611
+ speed_knots=float(v_turn),
1612
+ method=method,
1613
+ heading_final_deg=float(h1),
1614
  )
1615
+ pts_turn, _, _ = simulate_path(RADAR_FIX_LAT, RADAR_FIX_LON, turn_legs, method)
1616
+ lat0, lon0 = float(pts_turn[-1][0]), float(pts_turn[-1][1])
1617
+
1618
+ # initialize beam after turn
1619
+ beam = [{
1620
+ "lat": lat0,
1621
+ "lon": lon0,
1622
+ "headings": [],
1623
+ "speeds": [],
1624
+ "h_prev": float(h1),
1625
+ "v_turn": float(v_turn),
1626
+ "turn_legs": turn_legs,
1627
+ "turn_info": turn_info,
1628
+ "max_abs": 0.0,
1629
+ "sum_abs": 0.0,
1630
+ "smooth_sum": 0.0,
1631
+ "seg_end_errs": [],
1632
+ }]
1633
+
1634
+ # walk segments with beam search
1635
+ for mins, end_lbl in zip(GUIDED_SOUTH_SEG_MIN, GUIDED_SOUTH_SEG_LABELS_END):
1636
+ dt_end = _dt_from_label(end_lbl)
1637
+ meas = _lookup_measured_bto(dt_end)
1638
+ next_beam = []
1639
+
1640
+ for st in beam:
1641
+ h_prev = float(st["h_prev"])
1642
+ # candidate headings around previous
1643
+ cand_h = []
1644
+ for d in VAR_HEADING_DELTA_SET:
1645
+ h_try = h_prev + float(d)
1646
+ h_try = _clamp(h_try, VAR_HEADING_MIN, VAR_HEADING_MAX)
1647
+ cand_h.append(h_try)
1648
+ # unique small set
1649
+ cand_h = sorted(list({round(x, 6): x for x in cand_h}.values()))
1650
+
1651
+ for h_try in cand_h:
1652
+ v_seg, err_seg, (lat1, lon1) = _fit_segment_speed(
1653
+ method=method,
1654
+ lat0=float(st["lat"]),
1655
+ lon0=float(st["lon"]),
1656
+ hdg=float(h_try),
1657
+ minutes=float(mins),
1658
+ dt_end=dt_end,
1659
+ alt_ft=alt_ft,
1660
+ vs_fpm=vs_fpm,
1661
+ start_dt=start_dt,
1662
+ descent_start_dt=descent_start_dt,
1663
+ vmin=SOUTH_SPEED_MIN_KT,
1664
+ vmax=SOUTH_SPEED_MAX_KT,
1665
+ )
1666
+
1667
+ # BTO-first objective: update errors only if measured exists
1668
+ abs_err = 0.0 if meas is None or (err_seg != err_seg) else abs(float(err_seg))
1669
+ smooth_pen = abs(float(h_try) - float(h_prev)) * VAR_SMOOTH_PENALTY_US_PER_DEG
1670
+
1671
+ max_abs = max(float(st["max_abs"]), float(abs_err))
1672
+ sum_abs = float(st["sum_abs"]) + float(abs_err)
1673
+ smooth_sum = float(st["smooth_sum"]) + float(smooth_pen)
1674
+
1675
+ next_beam.append({
1676
+ "lat": float(lat1),
1677
+ "lon": float(lon1),
1678
+ "headings": st["headings"] + [float(h_try)],
1679
+ "speeds": st["speeds"] + [float(v_seg)],
1680
+ "h_prev": float(h_try),
1681
+ "v_turn": float(st["v_turn"]),
1682
+ "turn_legs": st["turn_legs"],
1683
+ "turn_info": st["turn_info"],
1684
+ "max_abs": max_abs,
1685
+ "sum_abs": sum_abs,
1686
+ "smooth_sum": smooth_sum,
1687
+ "seg_end_errs": st["seg_end_errs"] + [float(err_seg)],
1688
+ })
1689
+
1690
+ # prune beam
1691
+ if not next_beam:
1692
+ beam = []
1693
+ break
1694
+
1695
+ next_beam.sort(key=lambda x: (x["max_abs"], x["sum_abs"], x["smooth_sum"]))
1696
+ beam = next_beam[:max(1, int(VAR_BEAM_WIDTH))]
1697
+
1698
+ if not beam:
1699
  continue
1700
+
1701
+ # Evaluate final beam states using enforced timestamps (full recompute for safety)
1702
+ for st in beam:
1703
+ south_legs = _build_south_legs_variable(st["headings"], st["speeds"], GUIDED_SOUTH_SEG_MIN)
1704
+ full_legs = clip_legs_to_end_dt(st["turn_legs"] + south_legs, start_dt, end_dt)
1705
+
1706
+ bto_errs = _compute_bto_errors_at_enforced(method, start_dt, alt_ft, vs_fpm, descent_start_dt, full_legs, enforced_dts)
1707
+ if not bto_errs:
1708
+ continue
1709
+
1710
+ max_abs, sum_abs, rmse = score_from_errors(bto_errs)
1711
+ key = (max_abs, sum_abs, st["smooth_sum"], rmse)
1712
+
1713
+ if best_global is None or key < best_global:
1714
+ best_global = key
1715
+ best_global_legs = full_legs
1716
+ best_global_info = {
1717
+ "v_turn": float(st["v_turn"]),
1718
+ "turn_info": st["turn_info"],
1719
+ "headings": [float(x) for x in st["headings"]],
1720
+ "speeds": [float(x) for x in st["speeds"]],
1721
+ "seg_end_errs": [float(x) for x in st["seg_end_errs"]],
1722
+ "enforced_labels": enforced_labels,
1723
+ "max_abs_err_us": float(max_abs),
1724
+ "sum_abs_err_us": float(sum_abs),
1725
+ "rmse_bto_us": float(rmse),
1726
+ }
1727
+
1728
+ if best_global_legs is None:
1729
+ return None, "Guided solve: Variable mode found no solution within bounds.", None
1730
+
1731
+ hs = best_global_info["headings"]
1732
+ vs = best_global_info["speeds"]
1733
+ es = best_global_info["seg_end_errs"]
1734
 
1735
  report_lines = [
1736
+ "**Guided mode: Variable heading & speed + computed turn (18:22→18:40), clipped to 00:20**",
1737
+ "- Solver mode: **BTO-first ring-fit (beam search)**",
1738
+ f"- Turn speed: **{best_global_info['v_turn']:.2f} kt** (range {TURN_SPEED_MIN_KT:.0f}–{TURN_SPEED_MAX_KT:.0f})",
1739
+ f"- Beam width: **{VAR_BEAM_WIDTH}** | Heading deltas: **{VAR_HEADING_DELTA_SET}°** | Smooth penalty: **{VAR_SMOOTH_PENALTY_US_PER_DEG:.2f} µs/deg**",
 
1740
  "",
1741
+ "**Per-segment heading & speed:**",
1742
+ f"- 18:40→19:40: **{hs[0]:.2f}°**, **{vs[0]:.2f} kt** (end err {es[0]:+.2f} µs)",
1743
+ f"- 19:40→20:40: **{hs[1]:.2f}°**, **{vs[1]:.2f} kt** (end err {es[1]:+.2f} µs)",
1744
+ f"- 20:40→21:40: **{hs[2]:.2f}°**, **{vs[2]:.2f} kt** (end err {es[2]:+.2f} µs)",
1745
+ f"- 21:40→22:40: **{hs[3]:.2f}°**, **{vs[3]:.2f} kt** (end err {es[3]:+.2f} µs)",
1746
+ f"- 22:40→00:10: **{hs[4]:.2f}°**, **{vs[4]:.2f} kt** (end err {es[4]:+.2f} µs)",
1747
+ f"- 00:10→00:20: **{hs[5]:.2f}°**, **{vs[5]:.2f} kt** (end err {es[5]:+.2f} µs)",
1748
  "",
1749
+ f"- Enforced BTO timestamps: **{', '.join(best_global_info['enforced_labels'])}**",
1750
+ f"- Max |BTO error| at enforced times: **{best_global_info['max_abs_err_us']:.2f} µs**",
1751
+ f"- RMSE_BTO at enforced times: **{best_global_info['rmse_bto_us']:.2f} µs**",
1752
  "",
1753
  "*If 18:40 has no measured BTO in your satLoc, it is automatically not enforced.*",
1754
  ]
1755
+ return best_global_legs, "\n".join(report_lines), start_dt
 
1756
 
1757
 
1758
  # ================= RUNNERS =================
 
1794
  return gr.update(), "", html_map, stats_md, bfo_fig, bfo_md, alt_fig, gr.update(headers=headers, value=rows)
1795
 
1796
 
1797
+ def run_guided(guided_mode, start_lat, start_lon, method, alt_ft, vs_fpm, descent_start_choice):
1798
+ if guided_mode == "Straight heading":
1799
+ legs, report_md, start_dt = solve_guided_straight(method, alt_ft, vs_fpm, descent_start_choice)
 
 
 
 
 
 
 
 
1800
  else:
1801
+ legs, report_md, start_dt = solve_guided_variable(method, alt_ft, vs_fpm, descent_start_choice)
 
 
 
 
1802
 
1803
  if legs is None:
1804
  fig = go.Figure()
 
1844
  f"- **Initial track (fixed):** {RADAR_INITIAL_TRACK_DEG:.1f}°\n"
1845
  f"- **Final resting point:** **00:20:00Z** (trajectory clipped exactly at this time)\n\n"
1846
  f"**Turn speed bounds:** {TURN_SPEED_MIN_KT:.0f}–{TURN_SPEED_MAX_KT:.0f} kt\n"
1847
+ f"**Southern speed bounds:** {SOUTH_SPEED_MIN_KT:.0f}–{SOUTH_SPEED_MAX_KT:.0f} kt\n"
1848
+ f"**Straight mode:** Top-K polish (top-{TOPK_POLISH})\n"
1849
+ f"**Variable mode:** beam search (width {VAR_BEAM_WIDTH}) for speed"
1850
+ )
1851
+
1852
+ guided_mode = gr.Dropdown(
1853
+ ["Straight heading", "Variable heading & speed"],
1854
+ value="Straight heading",
1855
+ label="Guided scenario",
1856
  )
1857
 
1858
  with gr.Row():
 
1871
  label="When does descent start? (VS applies after this time)",
1872
  )
1873
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1874
  with gr.Row():
1875
  guided_btn = gr.Button("Generate: Guided scenario")
1876
  manual_btn = gr.Button("Simulate: Manual legs")
 
1917
 
1918
  guided_btn.click(
1919
  run_guided,
1920
+ inputs=[guided_mode, start_lat, start_lon, method, alt_ft, vs_fpm, descent_start_choice],
 
 
 
 
1921
  outputs=[legs_df, status_md, out_map, stats, bfo_plot, bfo_score, alt_plot, validation_df],
1922
  )
1923