From 76bf773a2374bd227191886abf65c9ac4c06ce27 Mon Sep 17 00:00:00 2001 From: Bryan Harter <41062454+bryan-harter@users.noreply.github.com> Date: Wed, 11 Nov 2026 18:34:18 -0700 Subject: [PATCH] Handle CoDICE SCI-LUT transitions across ancillary files and energy grids --- .../codice/ancillary.rst | 13 +++ imap_processing/codice/codice_l1a.py | 10 ++- imap_processing/codice/codice_l1b.py | 39 ++++----- imap_processing/codice/codice_l2.py | 31 +++---- imap_processing/codice/utils.py | 61 ++++++++++++-- .../tests/codice/test_process_by_table_id.py | 77 ++++++++++++++++- .../tests/codice/test_transition_energy.py | 84 +++++++++++++++++++ 7 files changed, 269 insertions(+), 46 deletions(-) create mode 100644 imap_processing/tests/codice/test_transition_energy.py diff --git a/docs/source/algorithm-code-documentation/codice/ancillary.rst b/docs/source/algorithm-code-documentation/codice/ancillary.rst index d276feff52..53fd85a446 100644 --- a/docs/source/algorithm-code-documentation/codice/ancillary.rst +++ b/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 diff --git a/imap_processing/codice/codice_l1a.py b/imap_processing/codice/codice_l1a.py index fba72d2886..35d04f1c77 100644 --- a/imap_processing/codice/codice_l1a.py +++ b/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") diff --git a/imap_processing/codice/codice_l1b.py b/imap_processing/codice/codice_l1b.py index ab2a2a883c..09e6358078 100644 --- a/imap_processing/codice/codice_l1b.py +++ b/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 [ diff --git a/imap_processing/codice/codice_l2.py b/imap_processing/codice/codice_l2.py index 4bb79c9f3e..033f21cff9 100644 --- a/imap_processing/codice/codice_l2.py +++ b/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") ) diff --git a/imap_processing/codice/utils.py b/imap_processing/codice/utils.py index 52527c0182..ecbbb433a9 100644 --- a/imap_processing/codice/utils.py +++ b/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( diff --git a/imap_processing/tests/codice/test_process_by_table_id.py b/imap_processing/tests/codice/test_process_by_table_id.py index 3fe28d21da..81f3d524c7 100644 --- a/imap_processing/tests/codice/test_process_by_table_id.py +++ b/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() diff --git a/imap_processing/tests/codice/test_transition_energy.py b/imap_processing/tests/codice/test_transition_energy.py new file mode 100644 index 0000000000..c7b20e5759 --- /dev/null +++ b/imap_processing/tests/codice/test_transition_energy.py @@ -0,0 +1,84 @@ +"""Verify transition-day energy support retains its time dependence.""" + +from unittest.mock import Mock + +import numpy as np +import pandas as pd +import pytest +import xarray as xr + +from imap_processing.codice import codice_l2, constants +from imap_processing.codice.codice_l1b import convert_to_rates + + +@pytest.mark.parametrize("varying", [False, True]) +def test_l1b_energy_matches_each_epoch(varying): + """L1B uses each epoch's voltage table, keeping normal days one-dimensional.""" + data = { + name: (("epoch", "esa_step"), np.ones((2, 2))) + for name in constants.LO_COUNTERS_SINGLES_VARIABLE_NAMES + } + for name in [ + "nso_half_spin", + "sw_bias_gain_mode", + "st_bias_gain_mode", + "spin_period", + "nso_esa_step", + "nso_spin_sector", + ]: + data[name] = ("epoch", [1, 1]) + data["acquisition_time_per_esa_step"] = (("epoch", "esa_step"), np.ones((2, 2))) + dataset = xr.Dataset(data, coords={"epoch": [10, 20], "esa_step": [0, 1]}) + dataset["k_factor"] = xr.DataArray([2.0], dims="k_factor") + voltage = [[100.0, 200.0], [300.0, 400.0]] if varying else [100.0, 200.0] + dims = ("epoch", "esa_step") if varying else ("esa_step",) + dataset["voltage_table"] = xr.DataArray(voltage, dims=dims) + result = convert_to_rates(dataset, "lo-counters-singles") + assert result["energy_per_charge"].dims == dims + np.testing.assert_allclose(result["energy_per_charge"], np.asarray(voltage) * 0.002) + assert result["energy_per_charge_label"].dims == dims + if varying: + assert result["energy_per_charge"].attrs["DEPEND_0"] == "epoch" + assert result["energy_per_charge_label"].values.tolist() == [ + ["0.200", "0.400"], + ["0.600", "0.800"], + ] + + +@pytest.mark.parametrize("varying", [False, True]) +def test_l2_preserves_time_dependent_energy(monkeypatch, varying): + """The L2 entry point must not swap an axis to a two-dimensional coordinate.""" + + dims = ("epoch", "esa_step") if varying else ("esa_step",) + values = np.array([[0.2, 0.4], [0.6, 0.8]]) if varying else np.array([0.2, 0.4]) + dataset = xr.Dataset( + { + "energy_per_charge": (dims, values, {"DEPEND_1": "esa_step"}), + "energy_per_charge_label": (dims, np.char.mod("%.3f", values)), + "esa_step_label": ("esa_step", ["0", "1"]), + "counts": (("epoch", "esa_step"), np.ones((2, 2))), + }, + coords={"epoch": [10, 20], "esa_step": [0, 1]}, + ) + if varying: + dataset["energy_per_charge"].attrs["DEPEND_0"] = "epoch" + dependencies = Mock() + dependencies.get_file_paths.return_value = [ + "imap_codice_l1b_lo-sw-species_20251218_v001.0001.cdf" + ] + monkeypatch.setattr(codice_l2, "load_cdf", lambda _: dataset) + monkeypatch.setattr(codice_l2, "get_geometric_factor_lut", lambda _: None) + monkeypatch.setattr(codice_l2, "compute_geometric_factors", lambda *_: None) + monkeypatch.setattr( + codice_l2, "get_efficiency_lut", lambda _: pd.DataFrame({"product": ["sw"]}) + ) + # Isolate coordinate handling; intensity math already uses named dimensions. + monkeypatch.setattr(codice_l2, "process_lo_species_intensity", lambda ds, *_: ds) + result = codice_l2.process_codice_l2("lo-sw-species", dependencies) + np.testing.assert_array_equal(result["energy_per_charge"], values) + if varying: + assert result["energy_per_charge"].dims == ("epoch", "esa_step") + assert result["counts"].dims == ("epoch", "esa_step") + assert result["energy_per_charge"].attrs["DEPEND_0"] == "epoch" + else: + assert result["counts"].dims == ("epoch", "energy_per_charge")