More robust MATLAB verification test

scripts/generate_matlab_verification_schedule.py

@@ -0,0 +1,110 @@
+import argparse
+import csv
+import tomllib
+from datetime import datetime
+from pathlib import Path
+
+REPO_ROOT = Path(__file__).resolve().parent.parent
+DEFAULT_MISSIONS = REPO_ROOT / 'tests/data/verification/legacy/missions.toml'
+DEFAULT_OUTPUT = REPO_ROOT / 'tests/data/verification/legacy/matlab-schedule.csv'
+DEFAULT_AIRPORTS = REPO_ROOT / 'tests/data/airports/airports.csv'
+DEFAULT_AIRPORTS_OUTPUT = (
+    REPO_ROOT / 'tests/data/verification/legacy/matlab-airports.csv'
+)
+
+
+def compact_time(value: datetime) -> str:
+    return value.strftime('%H%M')
+
+
+def main() -> None:
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument('--missions', type=Path, default=DEFAULT_MISSIONS)
+    parser.add_argument('--output', type=Path, default=DEFAULT_OUTPUT)
+    parser.add_argument('--airports', type=Path, default=DEFAULT_AIRPORTS)
+    parser.add_argument('--airports-output', type=Path, default=DEFAULT_AIRPORTS_OUTPUT)
+    args = parser.parse_args()
+
+    with args.missions.open('rb') as fp:
+        flights = tomllib.load(fp)['flight']
+
+    fieldnames = [
+        'depapt',
+        'arrapt',
+        'deptim',
+        'arrtim',
+        'days',
+        'inpacft',
+        'seats',
+        'efffrom',
+        'effto',
+        'NFlts',
+    ]
+    with args.output.open('w', newline='') as fp:
+        writer = csv.DictWriter(fp, fieldnames=fieldnames, lineterminator='\n')
+        writer.writeheader()
+        for flight in flights:
+            departure = datetime.fromisoformat(flight['departure'])
+            arrival = datetime.fromisoformat(flight['arrival'])
+            service_date = departure.strftime('%Y%m%d')
+            writer.writerow(
+                {
+                    'depapt': flight['origin'],
+                    'arrapt': flight['destination'],
+                    'deptim': compact_time(departure),
+                    'arrtim': compact_time(arrival),
+                    'days': departure.isoweekday(),
+                    'inpacft': flight['aircraft_type'],
+                    'seats': flight['seats'],
+                    'efffrom': service_date,
+                    'effto': service_date,
+                    'NFlts': 1,
+                }
+            )
+
+    airport_codes = {
+        code for flight in flights for code in (flight['origin'], flight['destination'])
+    }
+    with args.airports.open(newline='', encoding='utf-8') as fp:
+        airports = {
+            row['iata_code']: row
+            for row in csv.DictReader(fp)
+            if row['iata_code'] in airport_codes
+        }
+    missing = airport_codes - airports.keys()
+    if missing:
+        parser.error(f'missing airport data for: {", ".join(sorted(missing))}')
+
+    # Match the format used by MATLAB AEIC textscan call in
+    # readAllAirportData.m.
+    airport_fields = [
+        'Code',
+        'ICAO Code',
+        'Country Code',
+        'Latitude',
+        'Longitude',
+        'Elev_ft',
+        'Nrunways',
+        'Longest_runway_ft',
+    ]
+    with args.airports_output.open('w', newline='', encoding='utf-8') as fp:
+        writer = csv.DictWriter(fp, fieldnames=airport_fields, lineterminator='\n')
+        writer.writeheader()
+        for code in sorted(airport_codes):
+            airport = airports[code]
+            writer.writerow(
+                {
+                    'Code': code,
+                    'ICAO Code': airport['ident'],
+                    'Country Code': airport['iso_country'],
+                    'Latitude': airport['latitude_deg'],
+                    'Longitude': airport['longitude_deg'],
+                    'Elev_ft': airport['elevation_ft'] or 0,
+                    'Nrunways': 1,
+                    'Longest_runway_ft': 0,
+                }
+            )
+
+
+if __name__ == '__main__':
+    main()

scripts/run_legacy_verification.py

@@ -17,6 +17,8 @@
 from AEIC.storage import access_recorder, track_file_accesses
 from AEIC.trajectories import Trajectory
 from AEIC.types import Fuel, Species, SpeciesValues
+from AEIC.units import NAUTICAL_MILES_TO_METERS
+from AEIC.utils import GEOD
 from AEIC.verification.legacy import LegacyTrajectory
 from AEIC.verification.metrics import (
     ComparisonMetrics,
@@ -26,13 +28,9 @@
 
 TRAJ_FIELDS = [
     'flight_time',
-    'ground_distance',
-    'latitude',
-    'longitude',
     'altitude',
     'fuel_flow',
     'aircraft_mass',
-    'azimuth',
     'true_airspeed',
     'rate_of_climb',
 ]
@@ -50,6 +48,36 @@
 COMPARISON_FIELDS = TRAJ_FIELDS + ['trajectory_indices']
 
 SKIP_FINAL_POINT_FIELDS = set(['true_airspeed'])
+MAPE_PCT_TOL = 0.25
+GROUND_DISTANCE_MAE_M_TOL = 5.0 * NAUTICAL_MILES_TO_METERS
+POSITION_ROUTE_PCT_TOL = 0.25
+AZIMUTH_MAE_DEG_TOL = 0.25
+
+
+def position_error_pct(
+    legacy_traj: Trajectory, new_traj: Trajectory, route_distance: float
+) -> float:
+    """Mean WGS84 point separation as a percentage of route distance."""
+    _, _, distance = GEOD.inv(
+        legacy_traj.longitude,
+        legacy_traj.latitude,
+        new_traj.longitude,
+        new_traj.latitude,
+    )
+    return float(np.mean(distance) / route_distance * 100.0)
+
+
+def circular_mae_deg(reference: np.ndarray, actual: np.ndarray) -> float:
+    """Mean absolute heading error with 0/360-degree wraparound."""
+    difference = (actual - reference + 180.0) % 360.0 - 180.0
+    return float(np.mean(np.abs(difference)))
+
+
+def ground_distance_mae_m(legacy_traj: Trajectory, new_traj: Trajectory) -> float:
+    """Mean absolute cumulative ground-distance error in meters."""
+    return float(
+        np.mean(np.abs(legacy_traj.ground_distance - new_traj.ground_distance))
+    )
 
 
 def metrics_page(
@@ -312,7 +340,9 @@ def main(report_file) -> None:
         fuel = Fuel.model_validate(tomllib.load(fp))
 
     # Create a single trajectory builder to fly all missions.
-    builder = tb.LegacyBuilder(options=tb.Options(iterate_mass=False))
+    builder = tb.LegacyBuilder(
+        options=tb.Options(iterate_mass=False, use_weather=False)
+    )
 
     failed = []
 
@@ -340,7 +370,27 @@ def main(report_file) -> None:
             )
 
             # Record any metrics that are outside tolerance.
-            bad_metrics = out_of_tolerance(metrics, mape_pct_tol=0.25)
+            bad_metrics = out_of_tolerance(metrics, mape_pct_tol=MAPE_PCT_TOL)
+            ground_distance_error = ground_distance_mae_m(legacy_traj, new_traj)
+            if ground_distance_error > GROUND_DISTANCE_MAE_M_TOL:
+                ground_distance_error_nm = (
+                    ground_distance_error / NAUTICAL_MILES_TO_METERS
+                )
+                bad_metrics.append(
+                    f'ground_distance ({ground_distance_error_nm:.4f} nmi MAE)'
+                )
+            position_error = position_error_pct(
+                legacy_traj, new_traj, mission.gc_distance
+            )
+            if position_error > POSITION_ROUTE_PCT_TOL:
+                bad_metrics.append(
+                    f'position ({position_error:.4f}% of route distance)'
+                )
+            azimuth_error = circular_mae_deg(
+                legacy_traj.azimuth[:-1], new_traj.azimuth[:-1]
+            )
+            if azimuth_error > AZIMUTH_MAE_DEG_TOL:
+                bad_metrics.append(f'azimuth ({azimuth_error:.4f} deg MAE)')
             if len(bad_metrics) > 0:
                 failed.append((mission.label, bad_metrics))
 

src/AEIC/trajectories/builders/adjustable_legacy.py

@@ -426,7 +426,11 @@ def _fly_level_change(
         )
         if len(altitudes) == 0:
             return
-        if (flight_phase == FlightPhase.CLIMB and altitudes[-1] < final_altitude) or (
+        # Snap floating-point step accumulation to the exact endpoint
+        # to avoid a near-zero extra segment.
+        if np.isclose(altitudes[-1], final_altitude, rtol=0.0, atol=1e-9):
+            altitudes[-1] = final_altitude
+        elif (flight_phase == FlightPhase.CLIMB and altitudes[-1] < final_altitude) or (
             flight_phase == FlightPhase.DESCENT and altitudes[-1] > final_altitude
         ):
             altitudes = np.append(altitudes, final_altitude)
@@ -457,7 +461,7 @@ def _fly_level_change(
             else:
                 track_vector = self.weather.get_track_vector(
                     time=self.mission.departure,
-                    gt_point=self.ground_track.location(pt.ground_distance),
+                    gt_point=self.ground_track.step(pt.ground_distance, 0.0),
                     altitude=start_altitude,
                     horizontal_airspeed=horizontal_airspeed,
                     track_azimuth=pt.azimuth,

src/AEIC/trajectories/builders/legacy.py

@@ -300,7 +300,11 @@ def _fly_level_change(
             if flight_phase == FlightPhase.CLIMB
             else -self.altitude_step,
         )
-        if (flight_phase == FlightPhase.CLIMB and altitudes[-1] < final_altitude) or (
+        # Snap floating-point step accumulation to the exact endpoint
+        # to avoid a near-zero extra segment.
+        if np.isclose(altitudes[-1], final_altitude, rtol=0.0, atol=1e-9):
+            altitudes[-1] = final_altitude
+        elif (flight_phase == FlightPhase.CLIMB and altitudes[-1] < final_altitude) or (
             flight_phase == FlightPhase.DESCENT and altitudes[-1] > final_altitude
         ):
             altitudes = np.append(altitudes, final_altitude)
@@ -331,7 +335,7 @@ def _fly_level_change(
             else:
                 track_vector = self.weather.get_track_vector(
                     time=self.mission.departure,
-                    gt_point=self.ground_track.location(pt.ground_distance),
+                    gt_point=self.ground_track.step(pt.ground_distance, 0.0),
                     altitude=start_altitude,
                     horizontal_airspeed=horizontal_airspeed,
                     track_azimuth=pt.azimuth,

src/AEIC/verification/legacy.py

@@ -76,15 +76,19 @@ def trajectory(self) -> Trajectory:
         retval.longitude = self.df.long.values
         retval.altitude = self.df.alt.values * FEET_TO_METERS
         retval.ground_distance = self.df.horDist.values * NAUTICAL_MILES_TO_METERS
-        # Correct weird jumps in MATLAB azimuth output.
+        # Correct the near-180-degree direction reversals that MATLAB emits
+        # after an altitude-change step overshoots the destination airport.
         tmp_azimuth = self.df.az.values.copy()
         for i in range(1, len(tmp_azimuth)):
-            if tmp_azimuth[i] - tmp_azimuth[i - 1] > 180:
+            difference = tmp_azimuth[i] - tmp_azimuth[i - 1]
+            if 90 < difference < 270:
                 tmp_azimuth[i:] -= 180
-            elif tmp_azimuth[i - 1] - tmp_azimuth[i] > 180:
+            elif -270 < difference < -90:
                 tmp_azimuth[i:] += 180
         retval.azimuth = tmp_azimuth
         retval.true_airspeed = self.df.TAS.values
+        retval.ground_speed = self.df.groundSpeed.values
+        retval.heading = self.df.heading.values
         retval.rate_of_climb = self.df.roc_fpm.values * FPM_TO_MPS
         retval.aircraft_mass = self.df.acMass.values
         retval.fuel_flow = self.df.fuelFlow.values / MINUTES_TO_SECONDS

tests/data/airports/airports.csv

@@ -802,3 +802,7 @@
 "27238","ZYHB","large_airport","Harbin Taiping International Airport","45.623402","126.25","457","AS","CN","CN-23","Harbin","yes","ZYHB","HRB","ZYHB","","","https://en.wikipedia.org/wiki/Harbin_Taiping_International_Airport",""
 "27242","ZYTL","large_airport","Dalian Zhoushuizi International Airport","38.965719","121.538477","107","AS","CN","CN-21","Dalian (Ganjingzi)","yes","ZYTL","DLC","ZYTL","","","https://en.wikipedia.org/wiki/Dalian_Zhoushuizi_International_Airport","Dalian Air Base"
 "27243","ZYTX","large_airport","Shenyang Taoxian International Airport","41.639801","123.483002","198","AS","CN","CN-21","Hunnan, Shenyang","yes","ZYTX","SHE","ZYTX","","","https://en.wikipedia.org/wiki/Shenyang_Taoxian_International_Airport",""
+"2599","ENTC","large_airport","Tromsø Airport","69.683296","18.9189","31","EU","NO","NO-55","Tromsø","yes","ENTC","TOS","ENTC","","http://www.avinor.no/en/airport/tromso","https://en.wikipedia.org/wiki/Troms%C3%B8_Airport","Langnes, Tromsøya, Tromso"
+"4959","NFFN","large_airport","Nadi International Airport","-17.761822","177.437843","59","OC","FJ","FJ-W","Nadi","yes","NFFN","NAN","NFFN","","","https://en.wikipedia.org/wiki/Nadi_International_Airport",""
+"5488","PKMJ","large_airport","Marshall Islands International Airport","7.06511","171.271656","6","OC","MH","MH-MAJ","Majuro Atoll","yes","PKMJ","MAJ","PKMJ","MAJ","","https://en.wikipedia.org/wiki/Marshall_Islands_International_Airport",""
+"6184","SLLP","large_airport","El Alto International Airport","-16.510272","-68.189416","13355","SA","BO","BO-L","La Paz / El Alto","yes","SLLP","LPB","SLLP","","https://naabol.gob.bo/aeropuerto-internacional-el-alto/","https://en.wikipedia.org/wiki/El_Alto_International_Airport",""