CoDICE: handle SCI-LUT transitions across files and energy grids

docs/source/algorithm-code-documentation/codice/ancillary.rst

@@ -71,6 +71,19 @@ Summary
 The SCI-LUT
 -----------
 
+Transition days can contain packets with more than one ``Table_ID``. L1A indexes
+all supplied ``l1a-sci-lut`` dependencies by table ID and routes each packet group
+to its matching source file. Identical entries repeated in cumulative files are
+accepted; conflicting definitions of one ID are rejected. Missing IDs fail before
+processing that APID. Both files remain in the CLI's dependency provenance.
+Dependency generation must supply the preceding delivery as well as a delivery
+that becomes effective during the day; selecting a single file by date cannot
+represent an intraday onboard change. When a voltage table differs between groups,
+L1A retains it as an epoch-dependent support variable; L1B preserves that shape
+when computing energy per charge. L2 retains the ESA-step dimension for a
+time-dependent energy grid rather than swapping it to a two-dimensional axis.
+
+
 **[DOC]** Sections 5.1, 7 and 8. Originates as ``26850.03-SCI-LUT-01.xls`` with
 tabs ``Plan``, ``ESA Sweep``, ``Lo Stepping``, ``Views``, ``Collapse_Lo``,
 ``Collapse_Hi``, ``Data Products - Lo`` and ``Data Products - Hi``. The

imap_processing/codice/codice_l1a.py

@@ -28,7 +28,7 @@
 from imap_processing.codice.codice_l1a_lo_priority import l1a_lo_priority
 from imap_processing.codice.codice_l1a_lo_species import l1a_lo_species
 from imap_processing.codice.constants import CODICEAPID
-from imap_processing.codice.utils import process_by_table_id
+from imap_processing.codice.utils import index_sci_lut_files, process_by_table_id
 from imap_processing.utils import packet_file_to_datasets
 
 logger = logging.getLogger(__name__)
@@ -69,11 +69,15 @@ def process_l1a(  # noqa: PLR0912
         xtce_file,
     )
     datasets = []
+    lut_files = None
     for apid in datasets_by_apid:
         if apid not in [CODICEAPID.COD_LO_PHA, CODICEAPID.COD_HI_PHA]:
             # Get LUT file. Direct events do not need LUT
-            lut_file = dependency.get_file_paths(descriptor="l1a-sci-lut")
-            lut_file = lut_file[0]
+            if lut_files is None:
+                lut_files = index_sci_lut_files(
+                    dependency.get_file_paths(descriptor="l1a-sci-lut")
+                )
+            lut_file = lut_files
 
         if apid == CODICEAPID.COD_LO_SW_SPECIES_COUNTS:
             logger.info("Processing Lo SW Species Counts")

imap_processing/codice/codice_l1b.py

@@ -59,33 +59,30 @@ def convert_to_rates(
         # Calculate energy_per_charge using voltage_table and k_factor
         # 1e3 is to convert eV to keV
         energy_per_charge = (
-            dataset["voltage_table"].values * dataset["k_factor"].values * 1e-3
+            dataset["voltage_table"] * dataset["k_factor"].squeeze(drop=True) * 1e-3
         )
-        dataset["energy_per_charge"] = xr.DataArray(
-            energy_per_charge,
-            dims=[
-                "esa_step",
-            ],
-            attrs=cdf_attrs.get_variable_attributes(
-                "energy_per_charge", check_schema=False
-            )
+        energy_attrs = (
+            cdf_attrs.get_variable_attributes("energy_per_charge", check_schema=False)
             if cdf_attrs is not None
-            else {},  # the codice ialirt pipeline calls this function but
-            # energy_per_charge is an intermediate variable so
-            # it doesn't need attributes
+            else {}
         )
-        dataset["energy_per_charge_label"] = xr.DataArray(
-            np.array([f"{value:.3f}" for value in energy_per_charge]),
-            dims=[
-                "esa_step",
-            ],
-            attrs=cdf_attrs.get_variable_attributes(
+        label_attrs = (
+            cdf_attrs.get_variable_attributes(
                 "energy_per_charge_label", check_schema=False
             )
             if cdf_attrs is not None
-            else {},  # the codice ialirt pipeline calls this function but
-            # energy_per_charge_label is an intermediate variable so
-            # it does not need attributes
+            else {}
+        )
+        if "epoch" in energy_per_charge.dims:
+            for attrs in (energy_attrs, label_attrs):
+                attrs.update(DEPEND_0="epoch", DEPEND_1="esa_step")
+        dataset["energy_per_charge"] = xr.DataArray(
+            energy_per_charge.data, dims=energy_per_charge.dims, attrs=energy_attrs
+        )
+        dataset["energy_per_charge_label"] = xr.DataArray(
+            np.char.mod("%.3f", energy_per_charge.data),
+            dims=energy_per_charge.dims,
+            attrs=label_attrs,
         )
 
     if descriptor in [

imap_processing/codice/codice_l2.py

@@ -1504,21 +1504,22 @@ def process_codice_l2(
             # resolves to physical keV/e values instead of the esa_step index.
             # esa_step/esa_step_label need to keep their own esa_step
             # dimension though, so put them back after the swap.
-            esa_step = l2_dataset["esa_step"]
-            esa_step_label = l2_dataset["esa_step_label"]
-            l2_dataset = l2_dataset.swap_dims({"esa_step": "energy_per_charge"})
-            l2_dataset["esa_step"] = xr.DataArray(
-                esa_step.data, dims=["esa_step"], attrs=esa_step.attrs
-            )
-            l2_dataset["esa_step_label"] = xr.DataArray(
-                esa_step_label.data, dims=["esa_step"], attrs=esa_step_label.attrs
-            )
-            # Drop the DEPEND_1 inherited from L1B (pointing at esa_step) so
-            # energy_per_charge doesn't itself depend on esa_step.
-            l2_dataset["energy_per_charge"].attrs.pop("DEPEND_1", None)
-            l2_dataset["energy_per_charge_label"].attrs["DEPEND_1"] = (
-                "energy_per_charge"
-            )
+            if "epoch" not in l2_dataset["energy_per_charge"].dims:
+                esa_step = l2_dataset["esa_step"]
+                esa_step_label = l2_dataset["esa_step_label"]
+                l2_dataset = l2_dataset.swap_dims({"esa_step": "energy_per_charge"})
+                l2_dataset["esa_step"] = xr.DataArray(
+                    esa_step.data, dims=["esa_step"], attrs=esa_step.attrs
+                )
+                l2_dataset["esa_step_label"] = xr.DataArray(
+                    esa_step_label.data, dims=["esa_step"], attrs=esa_step_label.attrs
+                )
+                # Drop the DEPEND_1 inherited from L1B (pointing at esa_step) so
+                # energy_per_charge doesn't itself depend on esa_step.
+                l2_dataset["energy_per_charge"].attrs.pop("DEPEND_1", None)
+                l2_dataset["energy_per_charge_label"].attrs["DEPEND_1"] = (
+                    "energy_per_charge"
+                )
             l2_dataset.attrs.update(
                 cdf_attrs.get_global_attributes("imap_codice_l2_lo-sw-species")
             )

imap_processing/codice/utils.py

@@ -90,6 +90,35 @@ def read_sci_lut(file_path: Path, table_id: str) -> dict:
     return sci_lut_data
 
 
+def index_sci_lut_files(file_paths: list[Path]) -> dict[str, Path]:
+    """Index supplied SCI-LUT files by table ID, rejecting conflicting entries.
+
+    Parameters
+    ----------
+    file_paths : list[pathlib.Path]
+        Ancillary files supplied as processing dependencies.
+
+    Returns
+    -------
+    dict[str, pathlib.Path]
+        Source file for each table ID. Identical repeated entries are allowed.
+    """
+    tables: dict[str, dict] = {}
+    sources: dict[str, Path] = {}
+    for path in file_paths:
+        for table_id, table in json.loads(path.read_text()).items():
+            if table_id == "description":
+                continue
+            if table_id in tables and tables[table_id] != table:
+                raise ValueError(
+                    f"Conflicting SCI-LUT definitions for table ID {table_id}: "
+                    f"{sources[table_id].name} and {path.name}."
+                )
+            tables[table_id] = table
+            sources.setdefault(table_id, path)
+    return sources
+
+
 def get_view_tab_info(json_data: dict, view_id: int, apid: int) -> dict:
     """
     Get the view table information for a specific view and APID.
@@ -156,7 +185,7 @@ def get_view_tab_obj(
 
 def process_by_table_id(
     unpacked_dataset: xr.Dataset,
-    lut_file: Path,
+    lut_file: Path | dict[str, Path],
     process_fn: Callable[..., xr.Dataset],
 ) -> xr.Dataset:
     """
@@ -173,8 +202,8 @@ def process_by_table_id(
     ----------
     unpacked_dataset : xarray.Dataset
         Full unpacked dataset from the L0 packet file.
-    lut_file : pathlib.Path
-        Path to the SCI-LUT JSON file passed through to *process_fn*.
+    lut_file : pathlib.Path or dict[str, pathlib.Path]
+        SCI-LUT file, or source files indexed by table ID for transition days.
     process_fn : Callable
         The private ``_process_xxx`` function to call for each table_id group.
         It must accept the signature
@@ -192,12 +221,21 @@ def process_by_table_id(
     plan_step = unpacked_dataset["plan_step"].values[0]
 
     unique_table_ids = np.unique(unpacked_dataset["table_id"].values)
+    if isinstance(lut_file, dict):
+        missing = sorted(
+            {str(table_id) for table_id in unique_table_ids} - lut_file.keys()
+        )
+        if missing:
+            raise ValueError(
+                "Supplied SCI-LUT files do not have data for table IDs: "
+                + ", ".join(missing)
+            )
     processed = [
         process_fn(
             unpacked_dataset.isel(
                 epoch=unpacked_dataset["table_id"].values == table_id
             ),
-            lut_file,
+            lut_file[str(table_id)] if isinstance(lut_file, dict) else lut_file,
             table_id,
             view_id,
             apid,
@@ -212,13 +250,24 @@ def process_by_table_id(
     # instead of expanding them along the epoch dimension.
     # Eg. voltage_table and k_factor are 1-D arrays that apply to all
     # epochs in the same way.
-    return xr.concat(
+    # Voltage tables can change with the onboard LUT. Only promote this support
+    # variable when it differs; keep identical tables one-dimensional.
+    varying = [
+        name
+        for name in ("voltage_table",)
+        if name in processed[0]
+        and any(not processed[0][name].equals(ds[name]) for ds in processed[1:])
+    ]
+    combined = xr.concat(
         processed,
         dim="epoch",
-        data_vars="minimal",
+        data_vars=varying or "minimal",
         coords="minimal",
         compat="equals",
     ).sortby("epoch")
+    for name in varying:
+        combined[name].attrs["DEPEND_0"] = "epoch"
+    return combined
 
 
 def get_collapse_pattern_shape(

imap_processing/tests/codice/test_process_by_table_id.py

@@ -1,11 +1,18 @@
 """Unit tests for CoDICE table-id grouping helpers."""
 
+import json
 from pathlib import Path
+from unittest.mock import Mock
 
 import numpy as np
+import pytest
 import xarray as xr
 
-from imap_processing.codice.utils import process_by_table_id
+from imap_processing.codice.utils import (
+    index_sci_lut_files,
+    process_by_table_id,
+    read_sci_lut,
+)
 
 
 def test_process_by_table_id_keeps_voltage_table_1d():
@@ -39,3 +46,71 @@ def _process_fn(group_ds, lut_file, table_id, view_id, apid, plan_id, plan_step)
     assert processed["epoch"].values.tolist() == [10, 20]
     assert processed["voltage_table"].dims == ("esa_step",)
     np.testing.assert_array_equal(processed["voltage_table"].values, voltage_table)
+
+
+def test_transition_uses_each_tables_source(tmp_path):
+    """A day's old and new packets must use their own ancillary definitions."""
+
+    old = tmp_path / "old.json"
+    new = tmp_path / "new.json"
+    old.write_text(json.dumps({"description": "old", "3972558773": {"factor": 2}}))
+    new.write_text(json.dumps({"description": "new", "3974360784": {"factor": 3}}))
+    sources = index_sci_lut_files([new, old])
+    ds = xr.Dataset(
+        {
+            "view_id": ("epoch", [0, 0, 0]),
+            "pkt_apid": ("epoch", [1156, 1156, 1156]),
+            "plan_id": ("epoch", [0, 0, 0]),
+            "plan_step": ("epoch", [0, 0, 0]),
+            "table_id": ("epoch", [3974360784, 3972558773, 3972558773]),
+            "sample": ("epoch", [10, 20, 30]),
+        },
+        coords={"epoch": [30, 10, 20]},
+    )
+
+    def process(group, lut, table_id, *_):
+        factor = read_sci_lut(lut, table_id)["factor"]
+        return xr.Dataset(
+            {
+                "sample": group["sample"] * factor,
+                "voltage_table": ("esa_step", [factor, factor * 2]),
+            },
+            coords={"epoch": group.epoch},
+        )
+
+    result = process_by_table_id(ds, sources, process)
+    assert result.epoch.values.tolist() == [10, 20, 30]
+    assert result["sample"].values.tolist() == [40, 60, 30]
+    assert result["voltage_table"].dims == ("epoch", "esa_step")
+    assert result["voltage_table"].values.tolist() == [[2, 4], [2, 4], [3, 6]]
+    assert result["voltage_table"].attrs["DEPEND_0"] == "epoch"
+
+
+def test_lut_index_identical_duplicates_and_conflicts(tmp_path):
+    """Cumulative files may repeat a definition, but conflicting copies fail."""
+
+    old, new = tmp_path / "old.json", tmp_path / "new.json"
+    old.write_text(json.dumps({"description": "old", "1": {"factor": 2}}))
+    new.write_text(json.dumps({"description": "new", "1": {"factor": 2}, "2": {}}))
+    assert index_sci_lut_files([old, new]) == {"1": old, "2": new}
+    new.write_text(json.dumps({"1": {"factor": 3}}))
+    with pytest.raises(
+        ValueError, match="Conflicting SCI-LUT definitions for table ID 1"
+    ):
+        index_sci_lut_files([old, new])
+
+
+def test_transition_missing_id_fails_before_processing():
+    """Report all missing IDs without partly processing the dataset."""
+
+    ds = xr.Dataset(
+        {
+            name: ("epoch", [1, 2])
+            for name in ["view_id", "pkt_apid", "plan_id", "plan_step", "table_id"]
+        },
+        coords={"epoch": [0, 1]},
+    )
+    process = Mock()
+    with pytest.raises(ValueError, match="table IDs: 2"):
+        process_by_table_id(ds, {"1": Path("old.json")}, process)
+    process.assert_not_called()