Alignment with legacy MATLAB code (AEIC v2.1) (#106)

* Fix time interpolation behavior

* Bring trajectory calculations into line with MATLAB v2.1 results

* Fix constant emissions with custom fuel file; rename files

* update verification data with correct edb data

* updated hc co tests

* remove PM from matlab verification

* Regenerate golden test data

---------

Co-authored-by: aditeya <aditeya@mit.edu>

Author: ian-ross

notebooks/test-cases.ipynb

@@ -564,6 +564,151 @@
    "source": [
     "print(NOx_speciation_from_matlab())"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d3f82ea5",
+   "metadata": {},
+   "source": [
+    "# 4. Test HCCO function"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c838c3f8",
+   "metadata": {},
+   "source": [
+    "The Boeing Fuel flow 2 method (from the Boeing Fuel Flow Method 2 (BFFM2) from DuBois, D. & Paynter, G. (2006). *Fuel Flow Method 2 for Estimating Aircraft Emissions*. SAE Technical Paper 2006-01-1987) calculates the emissions index of HC, CO from fuel flow, ambient temperatures and pressures, as well as reference fuel flow and EI points that it interpolates from.\n",
+    "\n",
+    "The below implementation recreates while also using the same matlab code in AEICv2 (https://github.mit.edu/LAE/PW_EEI/blob/main/AEIC/Model_Files/hccoEIsfunc.m) or (/home/aditeya/SM_Thesis/PW_EEI/AEIC/Model_Files/hccoEIsfunc.m)\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "3b6f6947",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "def EI_HCCO(\n",
+    "    ff_eval,\n",
+    "    x_EI,\n",
+    "    ff_cal,\n",
+    "    Tamb,\n",
+    "    Pamb,\n",
+    "):\n",
+    "    ff = np.asarray(ff_eval, dtype=float)\n",
+    "    x_EI = np.asarray(x_EI, dtype=float).reshape(-1)\n",
+    "    ff_ref = np.asarray(ff_cal, dtype=float).reshape(-1)\n",
+    "\n",
+    "    if x_EI.size != 4 or ff_ref.size != 4:\n",
+    "        raise ValueError(\"x_EI and ff_cal must be length-4 arrays.\")\n",
+    "\n",
+    "    with np.errstate(divide=\"ignore\", invalid=\"ignore\"):\n",
+    "        slope = (np.log10(x_EI[1]) - np.log10(x_EI[0])) / (\n",
+    "            np.log10(ff_ref[1]) - np.log10(ff_ref[0])\n",
+    "        )\n",
+    "        x2 = np.log10(ff_ref[0])\n",
+    "        x3 = np.log10(x_EI[0])\n",
+    "        horz = (np.log10(x_EI[2]) + np.log10(x_EI[3])) / 2.0\n",
+    "\n",
+    "        intercept = (\n",
+    "            2.0 * np.log10(ff_ref[0]) * slope\n",
+    "            + np.log10(x_EI[2])\n",
+    "            + np.log10(x_EI[3])\n",
+    "            - 2.0 * np.log10(x_EI[0])\n",
+    "        ) / (2.0 * slope)\n",
+    "\n",
+    "        if intercept > np.log10(ff_ref[2]):\n",
+    "            intercept = np.log10(ff_ref[2])\n",
+    "        elif (intercept < np.log10(ff_ref[1])) and (slope < 0):\n",
+    "            horz = np.log10(x_EI[1])\n",
+    "            intercept = np.log10(ff_ref[1])\n",
+    "        elif slope >= 0:\n",
+    "            slope = 0.0\n",
+    "            x2 = 0.0\n",
+    "            x3 = horz\n",
+    "            intercept = np.log10(ff_ref[1])\n",
+    "\n",
+    "        xEI = np.zeros_like(ff, dtype=float)\n",
+    "        logff = np.log10(ff)\n",
+    "\n",
+    "        mask_slant = (logff < intercept) & (ff > 0)\n",
+    "        if np.any(mask_slant):\n",
+    "            xEI[mask_slant] = 10.0 ** (slope * (logff[mask_slant] - x2) + x3)\n",
+    "\n",
+    "        mask_horz = (logff >= intercept)\n",
+    "        if np.any(mask_horz):\n",
+    "            xEI[mask_horz] = 10.0 ** horz\n",
+    "\n",
+    "        xEI = np.where(np.isnan(xEI), 0.0, xEI)\n",
+    "\n",
+    "        theta_amb = np.asarray(Tamb, dtype=float) / 288.15\n",
+    "        delta_amb = np.asarray(Pamb, dtype=float) / 101325.0\n",
+    "        xEI = xEI * (theta_amb ** 3.3) / (delta_amb ** 1.02)\n",
+    "\n",
+    "    return xEI"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "ee378dfc",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# FROM EDB excel\n",
+    "EDB_fuel_flows = np.array([0.4, 0.8, 1.2, 1.8])\n",
+    "EDB_EI_HC = np.array([1.54, 0.05, 0.02, 0.03])\n",
+    "\n",
+    "\n",
+    "# Temperatures [K] and pressures [Pa] along test trajectory (calculated from ISA atmosphere).\n",
+    "temperatures = np.array([288.15, 278.4, 249.15, 216.65, 229.65, 275.15])\n",
+    "pressures = np.array([101325.0, 84555.9940737564, 47181.0021852292, 22632.0400950078, 30742.4326120969, 79495.201934051])\n",
+    "\n",
+    "# Total fuel flows along test trajectory [kg/s].\n",
+    "full_fuel_flows = np.array([0.3, 0.35, 0.55, 0.65, 0.5, 0.32])\n",
+    "\n",
+    "# Per-engine fuel flows along test trajectory [kg/s].\n",
+    "number_of_engines = 2\n",
+    "fuel_flows = [ff / number_of_engines for ff in full_fuel_flows]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "id": "64cf0c4f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "HC_result = EI_HCCO(\n",
+    "    fuel_flows,\n",
+    "    EDB_EI_HC,\n",
+    "    EDB_fuel_flows,\n",
+    "    Tamb=temperatures,\n",
+    "    Pamb=pressures,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "id": "e6db5e00",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[196.73236113  98.54714941  13.25415652   7.73938389  25.11776059\n",
+      " 157.25298236]\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(HC_result)"
+   ]
   }
  ],
  "metadata": {

scripts/run_legacy_verification.py

@@ -11,7 +11,7 @@
 from matplotlib.backends.backend_pdf import PdfPages
 
 import AEIC.trajectories.builders as tb
-from AEIC.config import Config, config
+from AEIC.config import Config
 from AEIC.emissions import compute_emissions
 from AEIC.missions import Mission
 from AEIC.performance.models import PerformanceModel
@@ -25,19 +25,22 @@
 )
 
 TRAJ_FIELDS = [
-    'fuel_flow',
-    'aircraft_mass',
+    'flight_time',
     'ground_distance',
     'latitude',
     'longitude',
+    'altitude',
+    'fuel_flow',
+    'aircraft_mass',
     'azimuth',
     'true_airspeed',
-    'altitude',
     'rate_of_climb',
 ]
 
 COMPARISON_FIELDS = TRAJ_FIELDS + ['trajectory_indices']
 
+SKIP_FINAL_POINT_FIELDS = set(['true_airspeed'])
+
 
 def metrics_page(
     pdf: PdfPages,
@@ -64,14 +67,19 @@ def metrics_page(
     # Title: mission ID.
     fig.text(0.5, 0.9, mission_id, ha='center', va='center', fontsize=24, weight='bold')
 
-    def t(x, line, txt):
+    def t(x, line, txt, error=False):
+        kwargs = {}
+        if error:
+            kwargs['color'] = 'red'
+            kwargs['weight'] = 'bold'
         fig.text(
             lmargin + x,
             tmargin - line * line_height,
             txt,
             ha='left',
             va='center',
             fontsize=14,
+            **kwargs,
         )
 
     # Display flight times.
@@ -82,7 +90,7 @@ def t(x, line, txt):
     t(0.5, 0, 'Legacy point count:')
     t(0.5, 1, 'New point count:')
     t(0.76, 0, f'{len(legacy_traj)}')
-    t(0.76, 1, f'{len(new_traj)}')
+    t(0.76, 1, f'{len(new_traj)}', error=len(new_traj) != len(legacy_traj))
 
     col_width = 0.1
 
@@ -103,10 +111,24 @@ def cell(col, row, txt, bold=False, left=False):
 
     # Table row headers.
     for idx, field in enumerate(TRAJ_FIELDS):
-        cell(0, idx + 1, field, bold=True, left=True)
+        cell(
+            0,
+            idx + 1,
+            field + ('*' if field in SKIP_FINAL_POINT_FIELDS else ''),
+            bold=True,
+            left=True,
+        )
     for idx, sp in enumerate(species):
         cell(0, idx + len(TRAJ_FIELDS) + 1, f'EI_{sp.name}', bold=True, left=True)
 
+    cell(
+        0,
+        len(TRAJ_FIELDS) + len(species) + 2,
+        "* Omit final point from comparison "
+        "(MATLAB code doesn't fill in the last point)",
+        left=True,
+    )
+
     # Table values.
     for idx, field in enumerate(TRAJ_FIELDS):
         ms = metrics[field]
@@ -245,8 +267,9 @@ def run(report_file) -> None:
     # Set up paths to test data.
     data_dir = Path(__file__).parent.parent / 'tests/data/verification/legacy'
     legacy_dir = data_dir / 'matlab-output'
-    missions_file = data_dir / 'legacy_verf_missions.toml'
-    perf_path = data_dir / 'legacy_verification.toml'
+    missions_file = data_dir / 'missions.toml'
+    fuel_file = data_dir / 'fuel.toml'
+    perf_path = data_dir / 'performance-model.toml'
 
     # Set up AEIC configuration and set the AEIC_PATH to include the test data
     # directory.
@@ -258,7 +281,7 @@ def run(report_file) -> None:
     with open(missions_file, 'rb') as fp:
         mission_dict = tomllib.load(fp)
     missions = Mission.from_toml(mission_dict)
-    with open(config.emissions.fuel_file, 'rb') as fp:
+    with open(fuel_file, 'rb') as fp:
         fuel = Fuel.model_validate(tomllib.load(fp))
 
     # Create a single trajectory builder to fly all missions.
@@ -277,14 +300,17 @@ def run(report_file) -> None:
             new_traj = builder.fly(pm, mission)
             new_traj.add_fields(compute_emissions(pm, fuel, new_traj))
 
-            # For comparison, we need to interpolate the new AEIC trajectory onto
-            # the same time points as the legacy trajectory, since they may not
-            # match.
-            interp_traj = new_traj.interpolate_time(legacy_traj.flight_time)
+            # For comparison, we do *not* interpolate the new AEIC trajectory
+            # onto the same time points as the legacy trajectory. The match
+            # should be close enough that we can compare corresponding points
+            # along the trajectories. The number of points in the trajectories
+            # should match exactly.
 
             # Compute comparison metrics.
             # dict[str, ComparisonMetrics | SpeciesValues[ComparisonMetrics]]
-            metrics = legacy_traj.compare(interp_traj, COMPARISON_FIELDS)
+            metrics = legacy_traj.compare(
+                new_traj, COMPARISON_FIELDS, SKIP_FINAL_POINT_FIELDS
+            )
 
             # Record any metrics that are outside tolerance.
             bad_metrics = out_of_tolerance(metrics, rtol=0.05, atol=1.0e-3)

src/AEIC/emissions/ei/hcco.py

@@ -136,12 +136,12 @@ def EI_HCCO(
     #       xEI_acrp = xEI * [1 + (–52) * (ff_eval – ff_cal[0])]
     #    Then overwrite those points with xEI_acrp.
     # ----------------------------------------------------------------------------
-    ACRP_slope = -52.0
-    low_thrust_mask = ff_eval < ff_cal[ThrustMode.IDLE]
-    if np.any(low_thrust_mask):
-        delta_ff = ff_eval[low_thrust_mask] - ff_cal[ThrustMode.IDLE]
-        xEI_acrp = xEI_out[low_thrust_mask] * (1.0 + ACRP_slope * delta_ff)
-        xEI_out[low_thrust_mask] = xEI_acrp
+    # ACRP_slope = -52.0
+    # low_thrust_mask = ff_eval < ff_cal[ThrustMode.IDLE]
+    # if np.any(low_thrust_mask):
+    #     delta_ff = ff_eval[low_thrust_mask] - ff_cal[ThrustMode.IDLE]
+    #     xEI_acrp = xEI_out[low_thrust_mask] * (1.0 + ACRP_slope * delta_ff)
+    #     xEI_out[low_thrust_mask] = xEI_acrp
 
     # ----------------------------------------------------------------------------
     # 7. Cruise correction:

src/AEIC/storage/container.py

@@ -330,6 +330,8 @@ def make_point(self, idx: int | None = None) -> Container:
         if idx is not None:
             if idx < -self._size or idx >= self._size:
                 raise IndexError('point index out of range')
+            if idx < 0:
+                idx += self._size
             for name in pt._data_dictionary:
                 if name in self._data:
                     val = self._data[name]
@@ -344,25 +346,27 @@ def _expand_capacity(self) -> None:
                 if isinstance(self._data[name], np.ndarray):
                     self._data[name] = np.resize(self._data[name], (self._capacity,))
 
-    def append(self, point: Container | None = None, **kwargs) -> None:
+    def append(
+        self, point: Container | None = None, replace: bool = False, **kwargs
+    ) -> None:
         if not self._extensible:
             raise ValueError('cannot append to fixed-size Container')
         if point is not None and len(kwargs) != 0:
             raise ValueError('cannot specify both point and keyword arguments')
         if point is not None:
-            self._append_point(point)
+            self._append_point(point, replace)
         else:
-            self._append_from_dict(kwargs)
+            self._append_from_dict(kwargs, replace)
 
-    def _append_point(self, point: Container) -> None:
+    def _append_point(self, point: Container, replace: bool = False) -> None:
         # TODO: Cache single point fieldset.
         if point._data_dictionary != self._data_dictionary.single_point():
             raise ValueError('cannot append point with different data dictionary')
         if len(point) != 1:
             raise ValueError('can only append a single point')
-        self._append_from_dict(point._data)
+        self._append_from_dict(point._data, replace)
 
-    def _append_from_dict(self, data: dict[str, Any]) -> None:
+    def _append_from_dict(self, data: dict[str, Any], replace: bool = False) -> None:
         # Fields should match those in single point field set.
         field_set = self._data_dictionary.single_point()
         fields = set(field_set.keys())
@@ -375,13 +379,23 @@ def _append_from_dict(self, data: dict[str, Any]) -> None:
             if extra:
                 raise ValueError(f'extra fields in appended point: {extra}')
 
-        # Increase container capacity if needed.
-        if self._size == self._capacity:
-            self._expand_capacity()
+        if not replace:
+            # Increase container capacity if needed.
+            if self._size == self._capacity:
+                self._expand_capacity()
 
-        for name, value in data.items():
-            # Check that the type of the assigned value can be safely cast to
-            # the field type and cast and assign the value if OK.
-            self._data[name][self._size] = field_set[name].convert_in(value, name, 0)
+            for name, value in data.items():
+                # Check that the type of the assigned value can be safely cast to
+                # the field type and cast and assign the value if OK.
+                self._data[name][self._size] = field_set[name].convert_in(
+                    value, name, 0
+                )
 
-        self._size += 1
+            self._size += 1
+        else:
+            for name, value in data.items():
+                # Check that the type of the assigned value can be safely cast to
+                # the field type and cast and assign the value if OK.
+                self._data[name][self._size - 1] = field_set[name].convert_in(
+                    value, name, 0
+                )