From 9d10d67041bb3f592a72d0200f24f0e6e0a5d9ac Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 17 Jun 2026 10:28:08 +0000 Subject: [PATCH 1/5] phase A of issue #43 --- src/zagg/config.py | 184 ++++++++++++++++++++++++++++++++++++++ src/zagg/processing.py | 122 +++++++++++++++++++------ tests/test_config.py | 112 +++++++++++++++++++++++ tests/test_processing.py | 186 +++++++++++++++++++++++++++++++++++++++ 4 files changed, 578 insertions(+), 26 deletions(-) diff --git a/src/zagg/config.py b/src/zagg/config.py index 9acd02a5..f623cf80 100644 --- a/src/zagg/config.py +++ b/src/zagg/config.py @@ -21,6 +21,16 @@ class DataSourceDict(TypedDict): coordinates: dict[str, str] variables: dict[str, str] quality_filter: NotRequired[dict] + filters: NotRequired[list[dict]] + + +# Structured-predicate comparison operators (issue #43). ``in``/``not_in`` take a +# ``values`` list; the rest take a scalar ``value``. These are the only +# pushdown-eligible filter language; an ``expression`` filter is a base-level-only, +# aggregation-time escape hatch that forfeits pushdown. +_SCALAR_OPS = frozenset({"eq", "ne", "ge", "le", "lt", "gt"}) +_SET_OPS = frozenset({"in", "not_in"}) +FILTER_OPS = _SCALAR_OPS | _SET_OPS @dataclass @@ -169,6 +179,9 @@ def validate_config(config: PipelineConfig) -> None: if "start_date" not in temporal or "end_date" not in temporal: raise ValueError("bounds.temporal requires start_date and end_date") + # Validate the structured filter list (issue #43, Phase A) + _validate_filters(config.data_source) + ds_vars = set(config.data_source.get("variables", {}).keys()) agg_vars = config.aggregation.get("variables", {}) @@ -218,6 +231,142 @@ def validate_config(config: PipelineConfig) -> None: _validate_expression_columns(name, meta["expression"], ds_vars) +def get_filters(config: PipelineConfig) -> list[dict]: + """Return the ordered list of normalized data-source filters (issue #43). + + Two filter languages coexist: + + - **Structured predicates** ``{level?, dataset, column?, op, value|values, + keep?}`` are machine-inspectable and are the only kind eligible for read + pushdown (Phase C). ``op`` is one of :data:`FILTER_OPS`; ``in``/``not_in`` + take ``values`` (a list), the rest take a scalar ``value``. ``column`` is an + integer selector into an N-D flag array (e.g. ATL03 ``signal_conf_ph``). + ``keep`` (default ``True``) keeps matching rows; ``keep: false`` drops them. + - **Expression** filters ``{expression: ""}`` are a base-level-only, + aggregation-time escape hatch that forfeits pushdown (opaque to the planner). + + The flat ``quality_filter: {dataset, value}`` is sugar synthesizing one + base-level ``op: eq`` structured filter, so the ATL06 path is unchanged. An + explicit ``filters:`` list, when present, is used as-is (the flat + ``quality_filter`` is then ignored). + + Each returned filter carries a normalized ``level`` (``None`` means the base + level) and, for structured predicates, an explicit ``keep`` bool. + + Parameters + ---------- + config : PipelineConfig + + Returns + ------- + list[dict] + """ + return filters_from_data_source(config.data_source) + + +def filters_from_data_source(data_source: dict) -> list[dict]: + """Normalize the filter list from a raw ``data_source`` dict. + + Shared by :func:`get_filters` and the read path (which only holds the + ``data_source`` mapping). See :func:`get_filters` for the schema. + """ + explicit = data_source.get("filters") + if explicit is not None: + return [_normalize_filter(f) for f in explicit] + qf = data_source.get("quality_filter") + if qf is not None: + return [ + { + "level": None, + "dataset": qf["dataset"], + "column": None, + "op": "eq", + "value": qf["value"], + "keep": True, + } + ] + return [] + + +def _normalize_filter(f: dict) -> dict: + """Normalize one raw filter dict into canonical form (see :func:`get_filters`).""" + if "expression" in f: + return {"level": f.get("level"), "expression": f["expression"]} + op = f["op"] + out = { + "level": f.get("level"), + "dataset": f["dataset"], + "column": f.get("column"), + "op": op, + "keep": bool(f.get("keep", True)), + } + if op in _SET_OPS: + out["values"] = list(f["values"]) + else: + out["value"] = f["value"] + return out + + +def _validate_filters(data_source: dict) -> None: + """Validate the ``filters`` list (and that a flat ``quality_filter`` is sane). + + Raises ``ValueError`` on: unknown op, missing ``dataset``, ``in``/``not_in`` + without a list ``values``, scalar ops without ``value``, non-int ``column``, + a non-base-level ``expression`` filter, or wrong ``value`` type. ``column`` is + required for the N-D flag case but cannot be checked against array rank here + (no data); rank checks happen at read time. + """ + filters = data_source.get("filters") + if filters is None: + return + if not isinstance(filters, list): + raise ValueError("data_source.filters must be a list") + for i, f in enumerate(filters): + if not isinstance(f, dict): + raise ValueError(f"filter[{i}] must be a mapping") + if "expression" in f: + if "op" in f or "dataset" in f: + raise ValueError( + f"filter[{i}]: 'expression' filters take no 'op'/'dataset' " + "(base-level aggregation-time escape hatch, no pushdown)" + ) + if f.get("level") is not None: + raise ValueError( + f"filter[{i}]: 'expression' filters are base-level only " + "(level must be omitted)" + ) + if not isinstance(f["expression"], str): + raise ValueError(f"filter[{i}]: 'expression' must be a string") + continue + if "dataset" not in f: + raise ValueError(f"filter[{i}]: structured filter requires 'dataset'") + op = f.get("op") + if op not in FILTER_OPS: + raise ValueError( + f"filter[{i}]: unknown op {op!r} (allowed: {sorted(FILTER_OPS)})" + ) + col = f.get("column") + if col is not None and not isinstance(col, int): + raise ValueError( + f"filter[{i}]: 'column' must be an integer index (got {col!r})" + ) + if op in _SET_OPS: + if not isinstance(f.get("values"), list): + raise ValueError(f"filter[{i}]: op {op!r} requires a 'values' list") + for v in f["values"]: + if not isinstance(v, (int, float)): + raise ValueError( + f"filter[{i}]: 'values' must be numeric (got {v!r})" + ) + else: + if "value" not in f: + raise ValueError(f"filter[{i}]: op {op!r} requires a scalar 'value'") + if not isinstance(f["value"], (int, float)): + raise ValueError( + f"filter[{i}]: 'value' must be numeric (got {f['value']!r})" + ) + + def _is_numeric(s: str) -> bool: """Check if a string is a numeric literal.""" try: @@ -484,3 +633,38 @@ def evaluate_expression(expression: str, columns: dict[str, np.ndarray]) -> floa **columns, } return float(eval(expression, ns)) # noqa: S307 + + +def evaluate_filter_expression( + expression: str, columns: dict[str, np.ndarray] +) -> np.ndarray: + """Evaluate a boolean filter expression to a per-row mask (issue #43). + + Like :func:`evaluate_expression` but returns the raw boolean array rather than + a scalar float — the base-level ``expression`` filter escape hatch (e.g. + ``"(h_li > 0) & (s_li < 1)"``). Uses the same restricted namespace. + + Parameters + ---------- + expression : str + Python boolean expression over numpy and column variables. + columns : dict[str, np.ndarray] + Mapping of column names to arrays. + + Returns + ------- + numpy.ndarray + Boolean mask. + """ + ns = { + "__builtins__": {}, + "np": np, + "numpy": np, + "len": len, + "float": float, + "int": int, + "abs": abs, + "sum": sum, + **columns, + } + return np.asarray(eval(expression, ns), dtype=bool) # noqa: S307 diff --git a/src/zagg/processing.py b/src/zagg/processing.py index 49aec1d9..0247a27e 100644 --- a/src/zagg/processing.py +++ b/src/zagg/processing.py @@ -16,7 +16,14 @@ from zarr import config, open_array from zarr.abc.store import Store -from zagg.config import PipelineConfig, default_config, get_agg_fields, get_data_vars +from zagg.config import ( + PipelineConfig, + default_config, + evaluate_filter_expression, + filters_from_data_source, + get_agg_fields, + get_data_vars, +) from zagg.schema import ProcessingMetadata logger = logging.getLogger(__name__) @@ -486,6 +493,48 @@ def _kernel_aggregate( # -- end EXPERIMENTAL kernel path --------------------------------------------- +_COMPARE = { + "eq": np.equal, + "ne": np.not_equal, + "ge": np.greater_equal, + "le": np.less_equal, + "lt": np.less, + "gt": np.greater, +} + + +def _predicate_mask(arr: np.ndarray, f: dict) -> np.ndarray: + """Build a 1-D boolean keep-mask for one structured predicate (issue #43). + + ``f`` is a normalized structured filter (see :func:`zagg.config.get_filters`): + ``{op, column, value|values, keep}``. An integer ``column`` selects a column + from a 2-D flag array before comparing; it is required for N-D arrays and + rejected for 1-D arrays. ``keep: false`` inverts the result (drop matches). + """ + column = f.get("column") + if arr.ndim > 1: + if column is None: + raise ValueError( + f"filter on '{f['dataset']}': N-D array requires an integer 'column'" + ) + arr = arr[:, column] + elif column is not None: + raise ValueError( + f"filter on '{f['dataset']}': 'column' set but array is 1-D" + ) + + op = f["op"] + if op == "in": + mask = np.isin(arr, f["values"]) + elif op == "not_in": + mask = ~np.isin(arr, f["values"]) + else: + mask = _COMPARE[op](arr, f["value"]) + if not f.get("keep", True): + mask = ~mask + return mask + + def _read_group( h5obj, group: str, data_source: dict, shard_key: int, grid, arrow: bool = False ): @@ -497,7 +546,9 @@ def _read_group( """ coordinates = data_source["coordinates"] variables = data_source["variables"] - quality_filter = data_source.get("quality_filter") + filters = filters_from_data_source(data_source) + structured = [f for f in filters if "expression" not in f] + expressions = [f for f in filters if "expression" in f] # Resolve coordinate paths coord_paths = [path.format(group=group) for path in coordinates.values()] @@ -523,46 +574,65 @@ def _read_group( min_idx = int(indices[0]) max_idx = int(indices[-1]) + 1 - # Build hyperslice dataset list: variables + optional quality filter + # Build hyperslice dataset list: variables + any structured-filter flag arrays. + # Read each distinct path once; flag datasets may coincide with a variable. datasets = [] - for path_template in variables.values(): - path = path_template.format(group=group) - datasets.append({"dataset": path, "hyperslice": [(min_idx, max_idx)]}) - - if quality_filter is not None: - qf_path = quality_filter["dataset"].format(group=group) - datasets.append({"dataset": qf_path, "hyperslice": [(min_idx, max_idx)]}) + paths_seen = set() + var_paths = {col: tmpl.format(group=group) for col, tmpl in variables.items()} + for path in var_paths.values(): + if path not in paths_seen: + datasets.append({"dataset": path, "hyperslice": [(min_idx, max_idx)]}) + paths_seen.add(path) + filter_paths = {id(f): f["dataset"].format(group=group) for f in structured} + for path in filter_paths.values(): + if path not in paths_seen: + datasets.append({"dataset": path, "hyperslice": [(min_idx, max_idx)]}) + paths_seen.add(path) data = h5obj.readDatasets(datasets) # Apply spatial mask to sliced data mask_sliced = mask_spatial[min_idx:max_idx] - # Apply quality filter if configured - if quality_filter is not None: - qf_path = quality_filter["dataset"].format(group=group) - q_flag = data[qf_path][mask_sliced] - quality_mask = q_flag == quality_filter["value"] - if np.sum(quality_mask) == 0: - return None - else: - quality_mask = None + # Combine structured predicates as ANDed keep-masks (issue #43). A single + # base-level ``op: eq`` filter (the ATL06 quality_filter sugar) reproduces the + # former path byte-for-byte. + keep_mask = None + for f in structured: + flag = data[filter_paths[id(f)]][mask_sliced] + fmask = _predicate_mask(flag, f) + keep_mask = fmask if keep_mask is None else (keep_mask & fmask) + if keep_mask is not None and np.sum(keep_mask) == 0: + return None - # Build dataframe + # Build dataframe (variables sliced to spatial mask, then to the keep-mask) leaf_sliced = leaf_ids[min_idx:max_idx][mask_sliced] data_dict = {} - for col_name, path_template in variables.items(): - path = path_template.format(group=group) + for col_name, path in var_paths.items(): values = data[path][mask_sliced] - if quality_mask is not None: - values = values[quality_mask] + if keep_mask is not None: + values = values[keep_mask] data_dict[col_name] = values - if quality_mask is not None: - data_dict["leaf_id"] = leaf_sliced[quality_mask] + if keep_mask is not None: + data_dict["leaf_id"] = leaf_sliced[keep_mask] else: data_dict["leaf_id"] = leaf_sliced + # Base-level ``expression`` filters: aggregation-time escape hatch, evaluated + # over the already-read variable columns (forfeits pushdown, issue #43). + for f in expressions: + cols = {c: data_dict[c] for c in variables if c in data_dict} + emask = evaluate_filter_expression(f["expression"], cols) + if emask.shape != data_dict["leaf_id"].shape: + raise ValueError( + f"expression filter {f['expression']!r} must yield a per-row " + f"boolean mask (got shape {emask.shape})" + ) + if np.sum(emask) == 0: + return None + data_dict = {k: v[emask] for k, v in data_dict.items()} + if arrow: import pyarrow as pa diff --git a/tests/test_config.py b/tests/test_config.py index 3a5122e1..6ffca6fc 100644 --- a/tests/test_config.py +++ b/tests/test_config.py @@ -14,6 +14,7 @@ get_child_order, get_coords, get_data_vars, + get_filters, get_store_path, load_config, load_config_from_dict, @@ -228,6 +229,117 @@ def test_missing_required_section(self): validate_config(cfg) +# --------------------------------------------------------------------------- +# Structured filters (issue #43, Phase A) +# --------------------------------------------------------------------------- + +def _cfg_with_filters(filters=None, quality_filter=None): + """Minimal valid config with a custom data_source filter spec.""" + ds = {"variables": {"h_li": "/{group}/h_li"}} + if filters is not None: + ds["filters"] = filters + if quality_filter is not None: + ds["quality_filter"] = quality_filter + return PipelineConfig( + data_source=ds, + aggregation={"variables": { + "h_min": {"function": "min", "source": "h_li", "dtype": "float32"}, + }}, + output={"grid": {"type": "healpix", "parent_order": 6, "child_order": 12}}, + ) + + +class TestFilters: + def test_quality_filter_synthesizes_base_eq(self, atl06_config): + filters = get_filters(atl06_config) + assert len(filters) == 1 + f = filters[0] + assert f["op"] == "eq" + assert f["level"] is None + assert f["column"] is None + assert f["keep"] is True + assert f["value"] == 0 + assert f["dataset"].endswith("atl06_quality_summary") + + def test_no_filters_returns_empty(self): + cfg = _cfg_with_filters() + assert get_filters(cfg) == [] + + def test_explicit_filters_override_quality_filter(self): + cfg = _cfg_with_filters( + filters=[{"dataset": "/{group}/conf", "column": 0, "op": "ne", "value": 0}], + quality_filter={"dataset": "/{group}/qs", "value": 0}, + ) + filters = get_filters(cfg) + assert len(filters) == 1 + assert filters[0]["column"] == 0 + assert filters[0]["op"] == "ne" + + def test_normalize_set_op_keeps_values_list(self): + cfg = _cfg_with_filters( + filters=[{"dataset": "/d", "op": "in", "values": [2, 3, 4]}] + ) + f = get_filters(cfg)[0] + assert f["values"] == [2, 3, 4] + assert "value" not in f + + def test_normalize_keep_drop(self): + cfg = _cfg_with_filters( + filters=[{"dataset": "/d", "op": "eq", "value": 1, "keep": False}] + ) + assert get_filters(cfg)[0]["keep"] is False + + def test_expression_filter_normalized(self): + cfg = _cfg_with_filters(filters=[{"expression": "h_li > 0"}]) + f = get_filters(cfg)[0] + assert f["expression"] == "h_li > 0" + assert f["level"] is None + + def test_unknown_op_rejected(self): + cfg = _cfg_with_filters(filters=[{"dataset": "/d", "op": "between", "value": 1}]) + with pytest.raises(ValueError, match="unknown op"): + validate_config(cfg) + + def test_column_must_be_int(self): + cfg = _cfg_with_filters( + filters=[{"dataset": "/d", "column": "land", "op": "ne", "value": 0}] + ) + with pytest.raises(ValueError, match="must be an integer"): + validate_config(cfg) + + def test_set_op_requires_values_list(self): + cfg = _cfg_with_filters(filters=[{"dataset": "/d", "op": "in", "value": 3}]) + with pytest.raises(ValueError, match="requires a 'values' list"): + validate_config(cfg) + + def test_scalar_op_requires_value(self): + cfg = _cfg_with_filters(filters=[{"dataset": "/d", "op": "eq"}]) + with pytest.raises(ValueError, match="requires a scalar 'value'"): + validate_config(cfg) + + def test_bad_value_type_rejected(self): + cfg = _cfg_with_filters(filters=[{"dataset": "/d", "op": "eq", "value": "x"}]) + with pytest.raises(ValueError, match="must be numeric"): + validate_config(cfg) + + def test_missing_dataset_rejected(self): + cfg = _cfg_with_filters(filters=[{"op": "eq", "value": 0}]) + with pytest.raises(ValueError, match="requires 'dataset'"): + validate_config(cfg) + + def test_expression_with_level_rejected(self): + cfg = _cfg_with_filters(filters=[{"expression": "h_li > 0", "level": "segment"}]) + with pytest.raises(ValueError, match="base-level only"): + validate_config(cfg) + + def test_expression_with_op_rejected(self): + cfg = _cfg_with_filters( + filters=[{"expression": "h_li > 0", "op": "eq", "dataset": "/d"}] + ) + with pytest.raises(ValueError, match="take no 'op'"): + validate_config(cfg) + + # --------------------------------------------------------------------------- # Helper accessors # --------------------------------------------------------------------------- diff --git a/tests/test_processing.py b/tests/test_processing.py index 88fe881d..795156b4 100644 --- a/tests/test_processing.py +++ b/tests/test_processing.py @@ -13,6 +13,8 @@ _group_columns, _kernel_able, _kernel_aggregate, + _predicate_mask, + _read_group, calculate_cell_statistics, process_shard, write_dataframe_to_zarr, @@ -645,3 +647,187 @@ def test_group_template_substitution(self): ds = default_config().data_source path = ds["coordinates"]["latitude"].format(group="gt2r") assert path == "/gt2r/land_ice_segments/latitude" + + +# --------------------------------------------------------------------------- +# Structured filters in the read path (issue #43, Phase A) +# --------------------------------------------------------------------------- + +class _FakeH5: + """Stub h5coro object: ``readDatasets`` returns canned arrays by path. + + Honors the ``hyperslice`` bound so sliced reads mirror the real driver. + """ + + def __init__(self, arrays): + self._arrays = arrays + + def readDatasets(self, datasets): # noqa: N802 (mirror real h5coro API) + out = {} + for d in datasets: + if isinstance(d, str): + out[d] = self._arrays[d] + continue + path = d["dataset"] + arr = self._arrays[path] + hs = d.get("hyperslice") + if hs is not None: + lo, hi = hs[0] + arr = arr[lo:hi] + out[path] = arr + return out + + +class _ShardGrid: + """Grid stub: leaf id == row index; every row maps to ``shard_key`` 0, + so the spatial filter keeps all rows and the structured filters are + exercised in isolation.""" + + @staticmethod + def assign(lats, lons): + return np.arange(len(lats)) + + @staticmethod + def shards_of(leaf_ids): + return np.zeros(len(leaf_ids), dtype=int) + + +class TestPredicateMask: + def test_scalar_ops_1d(self): + arr = np.array([0, 1, 2, 3, 0]) + assert _predicate_mask( + arr, {"dataset": "/d", "op": "eq", "value": 0, "column": None} + ).tolist() == [True, False, False, False, True] + assert _predicate_mask( + arr, {"dataset": "/d", "op": "ge", "value": 2, "column": None} + ).tolist() == [False, False, True, True, False] + + def test_set_ops(self): + arr = np.array([2, 3, 4, 5]) + assert _predicate_mask( + arr, {"dataset": "/d", "op": "in", "values": [2, 4]} + ).tolist() == [True, False, True, False] + assert _predicate_mask( + arr, {"dataset": "/d", "op": "not_in", "values": [2, 4]} + ).tolist() == [False, True, False, True] + + def test_keep_false_inverts(self): + arr = np.array([0, 1, 0]) + assert _predicate_mask( + arr, {"dataset": "/d", "op": "eq", "value": 0, "keep": False} + ).tolist() == [False, True, False] + + def test_nd_column_slicing(self): + # 2-D flag array (5 rows x 3 surface-type columns) + arr = np.array([[0, 9, 9], [-2, 9, 9], [1, 9, 9], [-2, 9, 9], [3, 9, 9]]) + # signal_conf_ph-style: column 0, != -2 + mask = _predicate_mask(arr, {"dataset": "/d", "column": 0, "op": "ne", "value": -2}) + assert mask.tolist() == [True, False, True, False, True] + + def test_nd_requires_column(self): + arr = np.zeros((3, 2)) + with pytest.raises(ValueError, match="requires an integer 'column'"): + _predicate_mask(arr, {"dataset": "/d", "op": "eq", "value": 0}) + + def test_column_on_1d_rejected(self): + arr = np.zeros(3) + with pytest.raises(ValueError, match="array is 1-D"): + _predicate_mask(arr, {"dataset": "/d", "column": 0, "op": "eq", "value": 0}) + + +class TestReadGroupFilters: + def _data_source(self, **extra): + ds = { + "coordinates": {"latitude": "/lat", "longitude": "/lon"}, + "variables": {"h": "/h"}, + } + ds.update(extra) + return ds + + def test_quality_filter_eq_path(self): + h5 = _FakeH5({ + "/lat": np.array([1.0, 2.0, 3.0, 4.0]), + "/lon": np.array([1.0, 2.0, 3.0, 4.0]), + "/h": np.array([10.0, 20.0, 30.0, 40.0], dtype=np.float32), + "/qs": np.array([0, 1, 0, 1]), + }) + ds = self._data_source(quality_filter={"dataset": "/qs", "value": 0}) + df = _read_group(h5, "gt1l", ds, 0, _ShardGrid()) + assert df["h"].tolist() == [10.0, 30.0] + + def test_quality_filter_byte_identical_to_manual_eq(self): + # The synthesized base eq filter must reproduce the legacy mask exactly. + h = np.array([10.0, 20.0, 30.0, 40.0, 50.0], dtype=np.float32) + qs = np.array([0, 1, 0, 0, 1]) + h5 = _FakeH5({ + "/lat": np.arange(5.0), "/lon": np.arange(5.0), "/h": h, "/qs": qs, + }) + ds = self._data_source(quality_filter={"dataset": "/qs", "value": 0}) + df = _read_group(h5, "gt1l", ds, 0, _ShardGrid()) + expected = h[qs == 0] + assert df["h"].to_numpy().tobytes() == expected.tobytes() + + def test_2d_signal_conf_filter(self): + conf = np.array([[0], [-2], [4], [-2], [3]]) # column 0, surface type + h5 = _FakeH5({ + "/lat": np.arange(5.0), "/lon": np.arange(5.0), + "/h": np.arange(5.0, dtype=np.float32), "/conf": conf, + }) + ds = self._data_source( + filters=[{"dataset": "/conf", "column": 0, "op": "ne", "value": -2}] + ) + df = _read_group(h5, "gt1l", ds, 0, _ShardGrid()) + assert df["h"].tolist() == [0.0, 2.0, 4.0] + + def test_multiple_anded_filters(self): + h5 = _FakeH5({ + "/lat": np.arange(5.0), "/lon": np.arange(5.0), + "/h": np.arange(5.0, dtype=np.float32), + "/conf": np.array([[5], [5], [0], [5], [5]]), + "/pod": np.array([0, 0, 0, 1, 0]), + }) + ds = self._data_source(filters=[ + {"dataset": "/conf", "column": 0, "op": "ne", "value": 0}, + {"dataset": "/pod", "op": "eq", "value": 0}, + ]) + df = _read_group(h5, "gt1l", ds, 0, _ShardGrid()) + # row2 dropped by conf==0, row3 dropped by pod==1 + assert df["h"].tolist() == [0.0, 1.0, 4.0] + + def test_in_op_integer_column(self): + h5 = _FakeH5({ + "/lat": np.arange(5.0), "/lon": np.arange(5.0), + "/h": np.arange(5.0, dtype=np.float32), + "/conf": np.array([[2], [0], [3], [1], [4]]), + }) + ds = self._data_source( + filters=[{"dataset": "/conf", "column": 0, "op": "in", "values": [2, 3, 4]}] + ) + df = _read_group(h5, "gt1l", ds, 0, _ShardGrid()) + assert df["h"].tolist() == [0.0, 2.0, 4.0] + + def test_expression_filter_base_level(self): + h5 = _FakeH5({ + "/lat": np.arange(5.0), "/lon": np.arange(5.0), + "/h": np.array([-1.0, 2.0, -3.0, 4.0, 5.0], dtype=np.float32), + }) + ds = self._data_source(filters=[{"expression": "h > 0"}]) + df = _read_group(h5, "gt1l", ds, 0, _ShardGrid()) + assert df["h"].tolist() == [2.0, 4.0, 5.0] + + def test_no_filter_keeps_all(self): + h5 = _FakeH5({ + "/lat": np.arange(3.0), "/lon": np.arange(3.0), + "/h": np.arange(3.0, dtype=np.float32), + }) + df = _read_group(h5, "gt1l", self._data_source(), 0, _ShardGrid()) + assert df["h"].tolist() == [0.0, 1.0, 2.0] + + def test_all_filtered_returns_none(self): + h5 = _FakeH5({ + "/lat": np.arange(3.0), "/lon": np.arange(3.0), + "/h": np.arange(3.0, dtype=np.float32), + "/qs": np.array([1, 1, 1]), + }) + ds = self._data_source(quality_filter={"dataset": "/qs", "value": 0}) + assert _read_group(h5, "gt1l", ds, 0, _ShardGrid()) is None From 7968bfa433237e10662de260715919a808397c80 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 17 Jun 2026 22:09:25 +0000 Subject: [PATCH 2/5] fold phase A review findings --- src/zagg/config.py | 29 ++---- src/zagg/processing.py | 22 ++--- tests/test_config.py | 39 +++++--- tests/test_processing.py | 205 ++++++++++++++++++++++++++++----------- 4 files changed, 194 insertions(+), 101 deletions(-) diff --git a/src/zagg/config.py b/src/zagg/config.py index 09b10a99..63a44950 100644 --- a/src/zagg/config.py +++ b/src/zagg/config.py @@ -415,8 +415,7 @@ def _validate_filters(data_source: dict) -> None: ) if f.get("level") is not None: raise ValueError( - f"filter[{i}]: 'expression' filters are base-level only " - "(level must be omitted)" + f"filter[{i}]: 'expression' filters are base-level only (level must be omitted)" ) if not isinstance(f["expression"], str): raise ValueError(f"filter[{i}]: 'expression' must be a string") @@ -425,29 +424,21 @@ def _validate_filters(data_source: dict) -> None: raise ValueError(f"filter[{i}]: structured filter requires 'dataset'") op = f.get("op") if op not in FILTER_OPS: - raise ValueError( - f"filter[{i}]: unknown op {op!r} (allowed: {sorted(FILTER_OPS)})" - ) + raise ValueError(f"filter[{i}]: unknown op {op!r} (allowed: {sorted(FILTER_OPS)})") col = f.get("column") - if col is not None and not isinstance(col, int): - raise ValueError( - f"filter[{i}]: 'column' must be an integer index (got {col!r})" - ) + if col is not None and (not isinstance(col, int) or isinstance(col, bool)): + raise ValueError(f"filter[{i}]: 'column' must be an integer index (got {col!r})") if op in _SET_OPS: if not isinstance(f.get("values"), list): raise ValueError(f"filter[{i}]: op {op!r} requires a 'values' list") for v in f["values"]: - if not isinstance(v, (int, float)): - raise ValueError( - f"filter[{i}]: 'values' must be numeric (got {v!r})" - ) + if not isinstance(v, (int, float)) or isinstance(v, bool): + raise ValueError(f"filter[{i}]: 'values' must be numeric (got {v!r})") else: if "value" not in f: raise ValueError(f"filter[{i}]: op {op!r} requires a scalar 'value'") - if not isinstance(f["value"], (int, float)): - raise ValueError( - f"filter[{i}]: 'value' must be numeric (got {f['value']!r})" - ) + if not isinstance(f["value"], (int, float)) or isinstance(f["value"], bool): + raise ValueError(f"filter[{i}]: 'value' must be numeric (got {f['value']!r})") def _is_numeric(s: str) -> bool: @@ -811,9 +802,7 @@ def evaluate_expression(expression: str, columns: dict[str, np.ndarray]) -> floa return float(_eval_expression_raw(expression, columns)) -def evaluate_filter_expression( - expression: str, columns: dict[str, np.ndarray] -) -> np.ndarray: +def evaluate_filter_expression(expression: str, columns: dict[str, np.ndarray]) -> np.ndarray: """Evaluate a boolean filter expression to a per-row mask (issue #43). Like :func:`evaluate_expression` but returns the raw boolean array rather than diff --git a/src/zagg/processing.py b/src/zagg/processing.py index e0ab3f03..80e24c31 100644 --- a/src/zagg/processing.py +++ b/src/zagg/processing.py @@ -168,8 +168,7 @@ def _has_vector_fields(config: PipelineConfig) -> bool: :func:`_arrow_column`). """ return any( - get_output_signature(meta)["kind"] == "vector" - for meta in get_agg_fields(config).values() + get_output_signature(meta)["kind"] == "vector" for meta in get_agg_fields(config).values() ) @@ -680,14 +679,10 @@ def _predicate_mask(arr: np.ndarray, f: dict) -> np.ndarray: column = f.get("column") if arr.ndim > 1: if column is None: - raise ValueError( - f"filter on '{f['dataset']}': N-D array requires an integer 'column'" - ) + raise ValueError(f"filter on '{f['dataset']}': N-D array requires an integer 'column'") arr = arr[:, column] elif column is not None: - raise ValueError( - f"filter on '{f['dataset']}': 'column' set but array is 1-D" - ) + raise ValueError(f"filter on '{f['dataset']}': 'column' set but array is 1-D") op = f["op"] if op == "in": @@ -701,9 +696,7 @@ def _predicate_mask(arr: np.ndarray, f: dict) -> np.ndarray: return mask -def _read_group( - h5obj, group: str, data_source: dict, shard_key: int, grid, arrow: bool = False -): +def _read_group(h5obj, group: str, data_source: dict, shard_key: int, grid, arrow: bool = False): """Read and spatially filter one HDF5 group. Returns a ``pandas.DataFrame`` (default) or, when ``arrow=True``, a @@ -789,7 +782,12 @@ def _read_group( # over the already-read variable columns (forfeits pushdown, issue #43). for f in expressions: cols = {c: data_dict[c] for c in variables if c in data_dict} - emask = evaluate_filter_expression(f["expression"], cols) + try: + emask = evaluate_filter_expression(f["expression"], cols) + except NameError as e: + raise NameError( + f"expression filter {f['expression']!r} references an undefined name: {e}" + ) from e if emask.shape != data_dict["leaf_id"].shape: raise ValueError( f"expression filter {f['expression']!r} must yield a per-row " diff --git a/tests/test_config.py b/tests/test_config.py index ffe9f66b..ebb9d342 100644 --- a/tests/test_config.py +++ b/tests/test_config.py @@ -255,6 +255,7 @@ def test_missing_required_section(self): # Structured filters (issue #43, Phase A) # --------------------------------------------------------------------------- + def _cfg_with_filters(filters=None, quality_filter=None): """Minimal valid config with a custom data_source filter spec.""" ds = {"variables": {"h_li": "/{group}/h_li"}} @@ -264,9 +265,11 @@ def _cfg_with_filters(filters=None, quality_filter=None): ds["quality_filter"] = quality_filter return PipelineConfig( data_source=ds, - aggregation={"variables": { - "h_min": {"function": "min", "source": "h_li", "dtype": "float32"}, - }}, + aggregation={ + "variables": { + "h_min": {"function": "min", "source": "h_li", "dtype": "float32"}, + } + }, output={"grid": {"type": "healpix", "parent_order": 6, "child_order": 12}}, ) @@ -298,17 +301,13 @@ def test_explicit_filters_override_quality_filter(self): assert filters[0]["op"] == "ne" def test_normalize_set_op_keeps_values_list(self): - cfg = _cfg_with_filters( - filters=[{"dataset": "/d", "op": "in", "values": [2, 3, 4]}] - ) + cfg = _cfg_with_filters(filters=[{"dataset": "/d", "op": "in", "values": [2, 3, 4]}]) f = get_filters(cfg)[0] assert f["values"] == [2, 3, 4] assert "value" not in f def test_normalize_keep_drop(self): - cfg = _cfg_with_filters( - filters=[{"dataset": "/d", "op": "eq", "value": 1, "keep": False}] - ) + cfg = _cfg_with_filters(filters=[{"dataset": "/d", "op": "eq", "value": 1, "keep": False}]) assert get_filters(cfg)[0]["keep"] is False def test_expression_filter_normalized(self): @@ -355,12 +354,28 @@ def test_expression_with_level_rejected(self): validate_config(cfg) def test_expression_with_op_rejected(self): - cfg = _cfg_with_filters( - filters=[{"expression": "h_li > 0", "op": "eq", "dataset": "/d"}] - ) + cfg = _cfg_with_filters(filters=[{"expression": "h_li > 0", "op": "eq", "dataset": "/d"}]) with pytest.raises(ValueError, match="take no 'op'"): validate_config(cfg) + def test_bool_column_rejected(self): + # bool is a subclass of int; filter column: true must be rejected. + cfg = _cfg_with_filters(filters=[{"dataset": "/d", "column": True, "op": "eq", "value": 0}]) + with pytest.raises(ValueError, match="must be an integer"): + validate_config(cfg) + + def test_bool_value_rejected(self): + # bool is a subclass of int; filter value: true must be rejected. + cfg = _cfg_with_filters(filters=[{"dataset": "/d", "op": "eq", "value": True}]) + with pytest.raises(ValueError, match="must be numeric"): + validate_config(cfg) + + def test_bool_in_values_rejected(self): + # bool elements in a 'values' list must be rejected. + cfg = _cfg_with_filters(filters=[{"dataset": "/d", "op": "in", "values": [0, True]}]) + with pytest.raises(ValueError, match="must be numeric"): + validate_config(cfg) + # --------------------------------------------------------------------------- # Helper accessors diff --git a/tests/test_processing.py b/tests/test_processing.py index 0a207239..94cce0e7 100644 --- a/tests/test_processing.py +++ b/tests/test_processing.py @@ -651,7 +651,7 @@ def fake_read_group(*args, **kwargs): monkeypatch.setattr("zagg.processing._read_group", fake_read_group) monkeypatch.setattr("zagg.processing.h5coro.H5Coro", lambda *a, **k: object()) # Avoid resolving a real h5coro driver (s3driver import / creds plumbing). - monkeypatch.setattr("zagg.processing._make_url_rewriter", lambda driver: (lambda u: u)) + monkeypatch.setattr("zagg.processing._make_url_rewriter", lambda driver: lambda u: u) def test_kernel_branch_matches_default_path(self, monkeypatch): """process_shard(handoff="arrow-kernel") agrees with the default path on the @@ -784,7 +784,7 @@ def _vector_cfg(): "params": {"minlength": 3}, }, } - } + }, ) def test_has_vector_fields(self): @@ -921,9 +921,7 @@ def _vector_cfg(): } from zagg.config import PipelineConfig - return PipelineConfig( - data_source=cfg.data_source, aggregation=agg, output=cfg.output - ) + return PipelineConfig(data_source=cfg.data_source, aggregation=agg, output=cfg.output) def test_vector_leaf_to_zarr_to_read(self): pytest.importorskip("pyarrow") @@ -1000,18 +998,17 @@ class _OneChunkGrid: dtype="float32", fill_value=np.float32("nan"), ) - edges = pa.FixedSizeListArray.from_arrays( - pa.array(np.arange(8.0, dtype="float32")), 4 - ) + edges = pa.FixedSizeListArray.from_arrays(pa.array(np.arange(8.0, dtype="float32")), 4) table = pa.table({"edges": edges}) with pytest.raises(ValueError, match="one whole chunk"): - write_dataframe_to_zarr( - table, store, grid=_OneChunkGrid(), chunk_idx=(0,) - ) + write_dataframe_to_zarr(table, store, grid=_OneChunkGrid(), chunk_idx=(0,)) + + # --------------------------------------------------------------------------- # Structured filters in the read path (issue #43, Phase A) # --------------------------------------------------------------------------- + class _FakeH5: """Stub h5coro object: ``readDatasets`` returns canned arrays by path. @@ -1063,9 +1060,12 @@ def test_scalar_ops_1d(self): def test_set_ops(self): arr = np.array([2, 3, 4, 5]) - assert _predicate_mask( - arr, {"dataset": "/d", "op": "in", "values": [2, 4]} - ).tolist() == [True, False, True, False] + assert _predicate_mask(arr, {"dataset": "/d", "op": "in", "values": [2, 4]}).tolist() == [ + True, + False, + True, + False, + ] assert _predicate_mask( arr, {"dataset": "/d", "op": "not_in", "values": [2, 4]} ).tolist() == [False, True, False, True] @@ -1104,12 +1104,14 @@ def _data_source(self, **extra): return ds def test_quality_filter_eq_path(self): - h5 = _FakeH5({ - "/lat": np.array([1.0, 2.0, 3.0, 4.0]), - "/lon": np.array([1.0, 2.0, 3.0, 4.0]), - "/h": np.array([10.0, 20.0, 30.0, 40.0], dtype=np.float32), - "/qs": np.array([0, 1, 0, 1]), - }) + h5 = _FakeH5( + { + "/lat": np.array([1.0, 2.0, 3.0, 4.0]), + "/lon": np.array([1.0, 2.0, 3.0, 4.0]), + "/h": np.array([10.0, 20.0, 30.0, 40.0], dtype=np.float32), + "/qs": np.array([0, 1, 0, 1]), + } + ) ds = self._data_source(quality_filter={"dataset": "/qs", "value": 0}) df = _read_group(h5, "gt1l", ds, 0, _ShardGrid()) assert df["h"].tolist() == [10.0, 30.0] @@ -1118,9 +1120,14 @@ def test_quality_filter_byte_identical_to_manual_eq(self): # The synthesized base eq filter must reproduce the legacy mask exactly. h = np.array([10.0, 20.0, 30.0, 40.0, 50.0], dtype=np.float32) qs = np.array([0, 1, 0, 0, 1]) - h5 = _FakeH5({ - "/lat": np.arange(5.0), "/lon": np.arange(5.0), "/h": h, "/qs": qs, - }) + h5 = _FakeH5( + { + "/lat": np.arange(5.0), + "/lon": np.arange(5.0), + "/h": h, + "/qs": qs, + } + ) ds = self._data_source(quality_filter={"dataset": "/qs", "value": 0}) df = _read_group(h5, "gt1l", ds, 0, _ShardGrid()) expected = h[qs == 0] @@ -1128,37 +1135,47 @@ def test_quality_filter_byte_identical_to_manual_eq(self): def test_2d_signal_conf_filter(self): conf = np.array([[0], [-2], [4], [-2], [3]]) # column 0, surface type - h5 = _FakeH5({ - "/lat": np.arange(5.0), "/lon": np.arange(5.0), - "/h": np.arange(5.0, dtype=np.float32), "/conf": conf, - }) - ds = self._data_source( - filters=[{"dataset": "/conf", "column": 0, "op": "ne", "value": -2}] + h5 = _FakeH5( + { + "/lat": np.arange(5.0), + "/lon": np.arange(5.0), + "/h": np.arange(5.0, dtype=np.float32), + "/conf": conf, + } ) + ds = self._data_source(filters=[{"dataset": "/conf", "column": 0, "op": "ne", "value": -2}]) df = _read_group(h5, "gt1l", ds, 0, _ShardGrid()) assert df["h"].tolist() == [0.0, 2.0, 4.0] def test_multiple_anded_filters(self): - h5 = _FakeH5({ - "/lat": np.arange(5.0), "/lon": np.arange(5.0), - "/h": np.arange(5.0, dtype=np.float32), - "/conf": np.array([[5], [5], [0], [5], [5]]), - "/pod": np.array([0, 0, 0, 1, 0]), - }) - ds = self._data_source(filters=[ - {"dataset": "/conf", "column": 0, "op": "ne", "value": 0}, - {"dataset": "/pod", "op": "eq", "value": 0}, - ]) + h5 = _FakeH5( + { + "/lat": np.arange(5.0), + "/lon": np.arange(5.0), + "/h": np.arange(5.0, dtype=np.float32), + "/conf": np.array([[5], [5], [0], [5], [5]]), + "/pod": np.array([0, 0, 0, 1, 0]), + } + ) + ds = self._data_source( + filters=[ + {"dataset": "/conf", "column": 0, "op": "ne", "value": 0}, + {"dataset": "/pod", "op": "eq", "value": 0}, + ] + ) df = _read_group(h5, "gt1l", ds, 0, _ShardGrid()) # row2 dropped by conf==0, row3 dropped by pod==1 assert df["h"].tolist() == [0.0, 1.0, 4.0] def test_in_op_integer_column(self): - h5 = _FakeH5({ - "/lat": np.arange(5.0), "/lon": np.arange(5.0), - "/h": np.arange(5.0, dtype=np.float32), - "/conf": np.array([[2], [0], [3], [1], [4]]), - }) + h5 = _FakeH5( + { + "/lat": np.arange(5.0), + "/lon": np.arange(5.0), + "/h": np.arange(5.0, dtype=np.float32), + "/conf": np.array([[2], [0], [3], [1], [4]]), + } + ) ds = self._data_source( filters=[{"dataset": "/conf", "column": 0, "op": "in", "values": [2, 3, 4]}] ) @@ -1166,27 +1183,101 @@ def test_in_op_integer_column(self): assert df["h"].tolist() == [0.0, 2.0, 4.0] def test_expression_filter_base_level(self): - h5 = _FakeH5({ - "/lat": np.arange(5.0), "/lon": np.arange(5.0), - "/h": np.array([-1.0, 2.0, -3.0, 4.0, 5.0], dtype=np.float32), - }) + h5 = _FakeH5( + { + "/lat": np.arange(5.0), + "/lon": np.arange(5.0), + "/h": np.array([-1.0, 2.0, -3.0, 4.0, 5.0], dtype=np.float32), + } + ) ds = self._data_source(filters=[{"expression": "h > 0"}]) df = _read_group(h5, "gt1l", ds, 0, _ShardGrid()) assert df["h"].tolist() == [2.0, 4.0, 5.0] def test_no_filter_keeps_all(self): - h5 = _FakeH5({ - "/lat": np.arange(3.0), "/lon": np.arange(3.0), - "/h": np.arange(3.0, dtype=np.float32), - }) + h5 = _FakeH5( + { + "/lat": np.arange(3.0), + "/lon": np.arange(3.0), + "/h": np.arange(3.0, dtype=np.float32), + } + ) df = _read_group(h5, "gt1l", self._data_source(), 0, _ShardGrid()) assert df["h"].tolist() == [0.0, 1.0, 2.0] def test_all_filtered_returns_none(self): - h5 = _FakeH5({ - "/lat": np.arange(3.0), "/lon": np.arange(3.0), - "/h": np.arange(3.0, dtype=np.float32), - "/qs": np.array([1, 1, 1]), - }) + h5 = _FakeH5( + { + "/lat": np.arange(3.0), + "/lon": np.arange(3.0), + "/h": np.arange(3.0, dtype=np.float32), + "/qs": np.array([1, 1, 1]), + } + ) ds = self._data_source(quality_filter={"dataset": "/qs", "value": 0}) assert _read_group(h5, "gt1l", ds, 0, _ShardGrid()) is None + + def test_filter_dataset_coincides_with_variable_path(self): + # Filter dataset path == variable path exercises the dedup branch + # (path must appear exactly once in the h5coro read list). + h5 = _FakeH5( + { + "/lat": np.arange(5.0), + "/lon": np.arange(5.0), + "/h": np.array([0.0, 1.0, 2.0, 3.0, 4.0], dtype=np.float32), + } + ) + ds = self._data_source(filters=[{"dataset": "/h", "op": "ge", "value": 2.0}]) + df = _read_group(h5, "gt1l", ds, 0, _ShardGrid()) + assert df["h"].tolist() == [2.0, 3.0, 4.0] + + def test_expression_filter_after_structured(self): + # Expression filter ANDed after a structured predicate. + h5 = _FakeH5( + { + "/lat": np.arange(5.0), + "/lon": np.arange(5.0), + "/h": np.array([1.0, 2.0, 3.0, 4.0, 5.0], dtype=np.float32), + "/qs": np.array([0, 0, 1, 0, 0]), + } + ) + ds = self._data_source( + filters=[ + {"dataset": "/qs", "op": "eq", "value": 0}, + {"expression": "h > 2"}, + ] + ) + df = _read_group(h5, "gt1l", ds, 0, _ShardGrid()) + # structured drops row2 (qs==1); expression keeps h>2 from remainder + assert df["h"].tolist() == [4.0, 5.0] + + def test_two_sequential_expression_filters(self): + # Two sequential expression filters both applied. + h5 = _FakeH5( + { + "/lat": np.arange(6.0), + "/lon": np.arange(6.0), + "/h": np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0], dtype=np.float32), + } + ) + ds = self._data_source( + filters=[ + {"expression": "h > 2"}, + {"expression": "h < 6"}, + ] + ) + df = _read_group(h5, "gt1l", ds, 0, _ShardGrid()) + assert df["h"].tolist() == [3.0, 4.0, 5.0] + + def test_expression_filter_undefined_name_raises(self): + # Expression referencing an undefined name re-raises as NameError. + h5 = _FakeH5( + { + "/lat": np.arange(3.0), + "/lon": np.arange(3.0), + "/h": np.arange(3.0, dtype=np.float32), + } + ) + ds = self._data_source(filters=[{"expression": "undefined_col > 0"}]) + with pytest.raises(NameError, match="undefined name"): + _read_group(h5, "gt1l", ds, 0, _ShardGrid()) From 749cb109af10e7a467030e9e90cbdacce36ce48e Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 17 Jun 2026 22:13:57 +0000 Subject: [PATCH 3/5] phase B of issue #43 --- src/zagg/config.py | 137 +++++++++++++++++++++++++++++++ src/zagg/processing.py | 114 ++++++++++++++++++++++++-- tests/test_config.py | 162 +++++++++++++++++++++++++++++++++++++ tests/test_processing.py | 168 +++++++++++++++++++++++++++++++++++++++ 4 files changed, 574 insertions(+), 7 deletions(-) diff --git a/src/zagg/config.py b/src/zagg/config.py index 63a44950..5d46f0b0 100644 --- a/src/zagg/config.py +++ b/src/zagg/config.py @@ -13,6 +13,40 @@ import zagg.configs +class LinkDict(TypedDict): + """Per-level link to the next coarser level (issue #43, Phase B). + + A *link* describes a contiguous-range parent->child tiling: each parent segment + ``p`` covers base-rate indices ``[index_beg[p] - index_base, ...`` for + ``count[p]`` children. ``index_base`` shifts the raw ``index_beg`` values so + that Python 0-based indexing into the base array is straightforward. + + ``reference_index`` is a reserved slot for a future explicit-index-array variant + (non-contiguous children per parent); leave it ``None`` for the contiguous case. + """ + + to: str # key of the coarser level in ``levels`` + index_beg: str # HDF5 path for the per-parent start index array + count: str # HDF5 path for the per-parent child count array + index_base: NotRequired[int] # subtracted from index_beg values (default 0) + reference_index: NotRequired[str | None] # reserved; must be None + + +class LevelDict(TypedDict): + """One hierarchical level in a multi-rate HDF5 source (issue #43, Phase B). + + A source may have several rates (e.g. ATL03 ``photons`` and ``segments``). + Each level declares its own ``path``, ``coordinates``, and ``variables``, + plus an optional ``link`` to a coarser parent level. The flat single-level + form (no ``levels``/``base_level`` keys in ``data_source``) stays first-class. + """ + + path: str # HDF5 group path template (may contain ``{group}``) + coordinates: list[str] # coordinate dataset names within ``path`` + variables: list[str] # variable dataset names within ``path`` + link: NotRequired[LinkDict | None] + + class DataSourceDict(TypedDict): """Type hints for the ``data_source`` section of a pipeline config.""" @@ -22,6 +56,11 @@ class DataSourceDict(TypedDict): variables: dict[str, str] quality_filter: NotRequired[dict] filters: NotRequired[list[dict]] + # Hierarchical multi-level form (issue #43, Phase B). When present, the flat + # ``coordinates``/``variables`` keys are still accepted for the base level but + # ``levels`` + ``base_level`` take precedence for the read path. + levels: NotRequired[dict[str, LevelDict]] + base_level: NotRequired[str] # Structured-predicate comparison operators (issue #43). ``in``/``not_in`` take a @@ -178,6 +217,9 @@ def validate_config(config: PipelineConfig) -> None: # Validate the structured filter list (issue #43, Phase A) _validate_filters(config.data_source) + # Validate hierarchical multi-level form (issue #43, Phase B) + _validate_levels(config.data_source) + ds_vars = set(config.data_source.get("variables", {}).keys()) agg_vars = config.aggregation.get("variables", {}) @@ -441,6 +483,101 @@ def _validate_filters(data_source: dict) -> None: raise ValueError(f"filter[{i}]: 'value' must be numeric (got {f['value']!r})") +def _validate_levels(data_source: dict) -> None: + """Validate the hierarchical ``levels``/``base_level`` form (issue #43, Phase B). + + Rules: + - ``base_level`` must name a key in ``levels``. + - ``link.to`` in each level must name another key in ``levels``. + - ``link.index_base`` must be a non-negative int when present. + - ``link.reference_index`` must be ``None`` when present (reserved slot). + - Only ``base_level`` may omit ``link`` (it has no coarser parent). + - Flat single-level form (no ``levels`` key) is always valid. + """ + levels = data_source.get("levels") + if levels is None: + return + if not isinstance(levels, dict) or not levels: + raise ValueError("data_source.levels must be a non-empty mapping") + base_level = data_source.get("base_level") + if base_level is None: + raise ValueError("data_source.base_level is required when levels is present") + if base_level not in levels: + raise ValueError( + f"data_source.base_level {base_level!r} is not a key in levels " + f"(available: {sorted(levels)})" + ) + level_keys = set(levels) + for name, lvl in levels.items(): + if not isinstance(lvl, dict): + raise ValueError(f"levels.{name} must be a mapping") + if "path" not in lvl: + raise ValueError(f"levels.{name}: 'path' is required") + link = lvl.get("link") + if link is None: + if name != base_level: + raise ValueError( + f"levels.{name}: non-base levels must have a 'link' " + f"(only {base_level!r} may omit it)" + ) + continue + if not isinstance(link, dict): + raise ValueError(f"levels.{name}.link must be a mapping") + for field_name in ("to", "index_beg", "count"): + if field_name not in link: + raise ValueError(f"levels.{name}.link: '{field_name}' is required") + unknown = set(link) - {"to", "index_beg", "count", "index_base", "reference_index"} + if unknown: + raise ValueError( + f"levels.{name}.link: unknown fields {sorted(unknown)} " + f"(allowed: to, index_beg, count, index_base, reference_index)" + ) + if link["to"] not in level_keys: + raise ValueError( + f"levels.{name}.link.to {link['to']!r} is not a key in levels " + f"(available: {sorted(level_keys)})" + ) + index_base = link.get("index_base", 0) + if not isinstance(index_base, int) or isinstance(index_base, bool) or index_base < 0: + raise ValueError( + f"levels.{name}.link.index_base must be a non-negative int (got {index_base!r})" + ) + ref = link.get("reference_index") + if ref is not None: + raise ValueError( + f"levels.{name}.link.reference_index is reserved and must be null/omitted " + f"(explicit index-array variant not yet implemented)" + ) + + +def get_levels(config: "PipelineConfig") -> dict | None: + """Return the ``levels`` mapping from the data source, or ``None`` if flat. + + Parameters + ---------- + config : PipelineConfig + + Returns + ------- + dict or None + """ + return config.data_source.get("levels") + + +def get_base_level(config: "PipelineConfig") -> str | None: + """Return the ``base_level`` key from the data source, or ``None`` if flat. + + Parameters + ---------- + config : PipelineConfig + + Returns + ------- + str or None + """ + return config.data_source.get("base_level") + + def _is_numeric(s: str) -> bool: """Check if a string is a numeric literal.""" try: diff --git a/src/zagg/processing.py b/src/zagg/processing.py index 80e24c31..5d4e1c1a 100644 --- a/src/zagg/processing.py +++ b/src/zagg/processing.py @@ -668,6 +668,47 @@ def _kernel_aggregate( } +def _expand_mask_to_base( + coarse_mask: np.ndarray, + index_beg_arr: np.ndarray, + count_arr: np.ndarray, + index_base: int, + total_base_size: int, +) -> np.ndarray: + """Expand a coarse-rate boolean mask to a base-rate boolean mask (issue #43, Phase B). + + Each coarse parent ``p`` covers base-rate rows + ``index_beg_arr[p] - index_base, ..., index_beg_arr[p] - index_base + count_arr[p] - 1``. + The contiguity assumption: ranges do not overlap and together tile the full base array. + + Parameters + ---------- + coarse_mask : np.ndarray + 1-D boolean array of length ``n_parents``. + index_beg_arr : np.ndarray + Per-parent start index into the base array (before ``index_base`` shift). + count_arr : np.ndarray + Per-parent child count (number of base-rate rows this parent covers). + index_base : int + Subtracted from ``index_beg_arr`` to get 0-based base indices. + total_base_size : int + Length of the output base-rate array. + + Returns + ------- + np.ndarray + 1-D boolean array of length ``total_base_size``. + """ + out = np.zeros(total_base_size, dtype=bool) + for p, keep in enumerate(coarse_mask): + if not keep: + continue + beg = int(index_beg_arr[p]) - index_base + cnt = int(count_arr[p]) + out[beg : beg + cnt] = True + return out + + def _predicate_mask(arr: np.ndarray, f: dict) -> np.ndarray: """Build a 1-D boolean keep-mask for one structured predicate (issue #43). @@ -702,11 +743,33 @@ def _read_group(h5obj, group: str, data_source: dict, shard_key: int, grid, arro Returns a ``pandas.DataFrame`` (default) or, when ``arrow=True``, a ``pyarrow.Table`` carrying the identical columns. Returns ``None`` when the group has no observations in this shard. + + Supports two modes (issue #43, Phase B): + + *Flat* (no ``levels``/``base_level`` in ``data_source``): unchanged from Phase A — + all structured filters are applied directly to base-rate data. + + *Hierarchical* (``levels`` + ``base_level`` present): structured filters whose + normalized ``level`` key names a non-base level are applied at coarse rate, then + expanded to base-rate via the level's ``link`` arrays (``_expand_mask_to_base``). + Base-level structured filters and expression filters are unchanged. """ coordinates = data_source["coordinates"] variables = data_source["variables"] filters = filters_from_data_source(data_source) - structured = [f for f in filters if "expression" not in f] + base_level_key = data_source.get("base_level") + levels = data_source.get("levels") + # Partition filters: base-level structured, coarse-level structured, expressions. + base_structured = [ + f + for f in filters + if "expression" not in f and (f.get("level") is None or f.get("level") == base_level_key) + ] + coarse_structured = [ + f + for f in filters + if "expression" not in f and f.get("level") is not None and f.get("level") != base_level_key + ] expressions = [f for f in filters if "expression" in f] # Resolve coordinate paths @@ -733,7 +796,40 @@ def _read_group(h5obj, group: str, data_source: dict, shard_key: int, grid, arro min_idx = int(indices[0]) max_idx = int(indices[-1]) + 1 - # Build hyperslice dataset list: variables + any structured-filter flag arrays. + # --- Coarse-level filter expansion (Phase B) --- + # For each filter whose level is not the base level, read the coarse-rate + # flag array from the declared level path, build a coarse mask, then expand + # to base-rate via the level link arrays. AND the results into ``cross_mask``. + cross_mask: np.ndarray | None = None + if coarse_structured and levels is not None: + for f in coarse_structured: + level_key = f["level"] + lvl = levels[level_key] + flag_path = f["dataset"].format(group=group) + # Read the coarse flag array (full level, no hyperslice — we need all parents + # to align with link arrays which are also full-length). + coarse_data = h5obj.readDatasets([{"dataset": flag_path}]) + coarse_arr = coarse_data[flag_path] + coarse_fmask = _predicate_mask(coarse_arr, f) + # Read the link arrays from this level. + link = lvl["link"] + index_base = int(link.get("index_base", 0)) + ibeg_path = link["index_beg"].format(group=group) + cnt_path = link["count"].format(group=group) + link_data = h5obj.readDatasets( + [ + {"dataset": ibeg_path}, + {"dataset": cnt_path}, + ] + ) + ibeg_arr = link_data[ibeg_path] + cnt_arr = link_data[cnt_path] + expanded = _expand_mask_to_base(coarse_fmask, ibeg_arr, cnt_arr, index_base, len(lats)) + cross_mask = expanded if cross_mask is None else (cross_mask & expanded) + if cross_mask is not None and np.sum(cross_mask[min_idx:max_idx]) == 0: + return None + + # Build hyperslice dataset list: variables + any base-level structured-filter arrays. # Read each distinct path once; flag datasets may coincide with a variable. datasets = [] paths_seen = set() @@ -742,7 +838,7 @@ def _read_group(h5obj, group: str, data_source: dict, shard_key: int, grid, arro if path not in paths_seen: datasets.append({"dataset": path, "hyperslice": [(min_idx, max_idx)]}) paths_seen.add(path) - filter_paths = {id(f): f["dataset"].format(group=group) for f in structured} + filter_paths = {id(f): f["dataset"].format(group=group) for f in base_structured} for path in filter_paths.values(): if path not in paths_seen: datasets.append({"dataset": path, "hyperslice": [(min_idx, max_idx)]}) @@ -753,14 +849,18 @@ def _read_group(h5obj, group: str, data_source: dict, shard_key: int, grid, arro # Apply spatial mask to sliced data mask_sliced = mask_spatial[min_idx:max_idx] - # Combine structured predicates as ANDed keep-masks (issue #43). A single - # base-level ``op: eq`` filter (the ATL06 quality_filter sugar) reproduces the - # former path byte-for-byte. + # Combine base-level structured predicates as ANDed keep-masks (issue #43). keep_mask = None - for f in structured: + for f in base_structured: flag = data[filter_paths[id(f)]][mask_sliced] fmask = _predicate_mask(flag, f) keep_mask = fmask if keep_mask is None else (keep_mask & fmask) + + # AND in the cross-level expanded mask, aligned to the sliced window. + if cross_mask is not None: + cross_sliced = cross_mask[min_idx:max_idx][mask_sliced] + keep_mask = cross_sliced if keep_mask is None else (keep_mask & cross_sliced) + if keep_mask is not None and np.sum(keep_mask) == 0: return None diff --git a/tests/test_config.py b/tests/test_config.py index ebb9d342..25fcbd6e 100644 --- a/tests/test_config.py +++ b/tests/test_config.py @@ -11,10 +11,12 @@ default_config, evaluate_expression, get_agg_fields, + get_base_level, get_child_order, get_coords, get_data_vars, get_filters, + get_levels, get_output_signature, get_store_path, load_config, @@ -377,6 +379,166 @@ def test_bool_in_values_rejected(self): validate_config(cfg) +# --------------------------------------------------------------------------- +# Hierarchical levels and link validation (issue #43, Phase B) +# --------------------------------------------------------------------------- + + +def _minimal_two_level_ds(**overrides): + """Return a minimal two-level data_source dict with one segment->photon link.""" + ds = { + "reader": "h5coro", + "groups": ["gt1l"], + "coordinates": {"latitude": "/gt1l/ph_lat", "longitude": "/gt1l/ph_lon"}, + "variables": {"h": "/gt1l/h_ph"}, + "base_level": "photons", + "levels": { + "photons": { + "path": "/{group}/heights", + "coordinates": ["lat_ph", "lon_ph"], + "variables": ["h_ph"], + "link": None, + }, + "segments": { + "path": "/{group}/geolocation", + "coordinates": ["reference_photon_lat", "reference_photon_lon"], + "variables": ["signal_conf_ph"], + "link": { + "to": "photons", + "index_beg": "/{group}/geolocation/ph_index_beg", + "count": "/{group}/geolocation/segment_ph_cnt", + }, + }, + }, + } + ds.update(overrides) + return ds + + +def _cfg_with_levels(**overrides): + ds = _minimal_two_level_ds(**overrides) + return PipelineConfig( + data_source=ds, + aggregation={"variables": {"count": {"function": "len", "source": "h", "dtype": "int32"}}}, + output={"grid": {"type": "healpix", "parent_order": 6, "child_order": 12}}, + ) + + +class TestLevelsValidation: + def test_valid_two_level_config(self): + validate_config(_cfg_with_levels()) + + def test_flat_form_still_valid(self, atl06_config): + # Flat form (no levels/base_level) must still pass. + assert get_levels(atl06_config) is None + assert get_base_level(atl06_config) is None + validate_config(atl06_config) + + def test_get_levels_and_base_level(self): + cfg = _cfg_with_levels() + levels = get_levels(cfg) + assert levels is not None + assert "photons" in levels + assert "segments" in levels + assert get_base_level(cfg) == "photons" + + def test_base_level_must_name_a_key(self): + cfg = _cfg_with_levels(base_level="nonexistent") + with pytest.raises(ValueError, match="not a key in levels"): + validate_config(cfg) + + def test_link_to_must_name_a_key(self): + ds = _minimal_two_level_ds() + ds["levels"]["segments"]["link"]["to"] = "nonexistent" + cfg = _cfg_with_levels(**ds) + with pytest.raises(ValueError, match="not a key in levels"): + validate_config(cfg) + + def test_link_to_must_name_a_key2(self): + # Build config directly to avoid _cfg_with_levels merging issues + ds = _minimal_two_level_ds() + ds["levels"]["segments"]["link"]["to"] = "nonexistent" + cfg = PipelineConfig( + data_source=ds, + aggregation={"variables": {"c": {"function": "len", "dtype": "int32"}}}, + output={"grid": {"type": "healpix", "parent_order": 6, "child_order": 12}}, + ) + with pytest.raises(ValueError, match="not a key in levels"): + validate_config(cfg) + + def test_link_missing_required_field(self): + ds = _minimal_two_level_ds() + del ds["levels"]["segments"]["link"]["count"] + cfg = PipelineConfig( + data_source=ds, + aggregation={"variables": {"c": {"function": "len", "dtype": "int32"}}}, + output={"grid": {"type": "healpix", "parent_order": 6, "child_order": 12}}, + ) + with pytest.raises(ValueError, match="'count' is required"): + validate_config(cfg) + + def test_link_unknown_field_rejected(self): + ds = _minimal_two_level_ds() + ds["levels"]["segments"]["link"]["bogus"] = "x" + cfg = PipelineConfig( + data_source=ds, + aggregation={"variables": {"c": {"function": "len", "dtype": "int32"}}}, + output={"grid": {"type": "healpix", "parent_order": 6, "child_order": 12}}, + ) + with pytest.raises(ValueError, match="unknown fields"): + validate_config(cfg) + + def test_base_level_without_link_ok(self): + # base_level is the only level allowed to have link: None. + cfg = _cfg_with_levels() + validate_config(cfg) + assert cfg.data_source["levels"]["photons"]["link"] is None + + def test_non_base_level_without_link_rejected(self): + ds = _minimal_two_level_ds() + ds["levels"]["segments"]["link"] = None + cfg = PipelineConfig( + data_source=ds, + aggregation={"variables": {"c": {"function": "len", "dtype": "int32"}}}, + output={"grid": {"type": "healpix", "parent_order": 6, "child_order": 12}}, + ) + with pytest.raises(ValueError, match="must have a 'link'"): + validate_config(cfg) + + def test_levels_missing_base_level_key_rejected(self): + ds = _minimal_two_level_ds() + del ds["base_level"] + cfg = PipelineConfig( + data_source=ds, + aggregation={"variables": {"c": {"function": "len", "dtype": "int32"}}}, + output={"grid": {"type": "healpix", "parent_order": 6, "child_order": 12}}, + ) + with pytest.raises(ValueError, match="base_level is required"): + validate_config(cfg) + + def test_index_base_must_be_nonneg_int(self): + ds = _minimal_two_level_ds() + ds["levels"]["segments"]["link"]["index_base"] = -1 + cfg = PipelineConfig( + data_source=ds, + aggregation={"variables": {"c": {"function": "len", "dtype": "int32"}}}, + output={"grid": {"type": "healpix", "parent_order": 6, "child_order": 12}}, + ) + with pytest.raises(ValueError, match="non-negative int"): + validate_config(cfg) + + def test_reference_index_must_be_none(self): + ds = _minimal_two_level_ds() + ds["levels"]["segments"]["link"]["reference_index"] = "/some/path" + cfg = PipelineConfig( + data_source=ds, + aggregation={"variables": {"c": {"function": "len", "dtype": "int32"}}}, + output={"grid": {"type": "healpix", "parent_order": 6, "child_order": 12}}, + ) + with pytest.raises(ValueError, match="reserved"): + validate_config(cfg) + + # --------------------------------------------------------------------------- # Helper accessors # --------------------------------------------------------------------------- diff --git a/tests/test_processing.py b/tests/test_processing.py index 94cce0e7..b8e12818 100644 --- a/tests/test_processing.py +++ b/tests/test_processing.py @@ -12,6 +12,7 @@ _build_groups, _build_output, _concat_and_group, + _expand_mask_to_base, _group_columns, _has_vector_fields, _iter_carrier_columns, @@ -1281,3 +1282,170 @@ def test_expression_filter_undefined_name_raises(self): ds = self._data_source(filters=[{"expression": "undefined_col > 0"}]) with pytest.raises(NameError, match="undefined name"): _read_group(h5, "gt1l", ds, 0, _ShardGrid()) + + +# --------------------------------------------------------------------------- +# _expand_mask_to_base and cross-level filter path (issue #43, Phase B) +# --------------------------------------------------------------------------- + + +class TestExpandMaskToBase: + def test_single_parent_kept(self): + # 1 parent keeps base rows 0-2 (3 photons). + coarse = np.array([True]) + ibeg = np.array([0]) + cnt = np.array([3]) + out = _expand_mask_to_base(coarse, ibeg, cnt, index_base=0, total_base_size=3) + np.testing.assert_array_equal(out, [True, True, True]) + + def test_single_parent_dropped(self): + coarse = np.array([False]) + ibeg = np.array([0]) + cnt = np.array([3]) + out = _expand_mask_to_base(coarse, ibeg, cnt, index_base=0, total_base_size=3) + np.testing.assert_array_equal(out, [False, False, False]) + + def test_two_parents_alternating(self): + # parent 0 -> base rows 0-1 (kept); parent 1 -> base rows 2-4 (dropped). + coarse = np.array([True, False]) + ibeg = np.array([0, 2]) + cnt = np.array([2, 3]) + out = _expand_mask_to_base(coarse, ibeg, cnt, index_base=0, total_base_size=5) + np.testing.assert_array_equal(out, [True, True, False, False, False]) + + def test_index_base_shift(self): + # HDF5 1-based indexing: index_beg values start at 1. + coarse = np.array([False, True]) + ibeg = np.array([1, 4]) # 1-based + cnt = np.array([3, 2]) + out = _expand_mask_to_base(coarse, ibeg, cnt, index_base=1, total_base_size=5) + # parent1 covers base rows 3 and 4 (ibeg=4-1=3, cnt=2). + np.testing.assert_array_equal(out, [False, False, False, True, True]) + + def test_empty_coarse_mask(self): + coarse = np.array([False, False, False]) + ibeg = np.array([0, 2, 5]) + cnt = np.array([2, 3, 1]) + out = _expand_mask_to_base(coarse, ibeg, cnt, index_base=0, total_base_size=6) + assert not np.any(out) + + def test_full_coarse_mask(self): + coarse = np.array([True, True]) + ibeg = np.array([0, 3]) + cnt = np.array([3, 2]) + out = _expand_mask_to_base(coarse, ibeg, cnt, index_base=0, total_base_size=5) + assert np.all(out) + + +class TestReadGroupCrossLevel: + """Phase B: cross-level filters expand coarse verdicts to base-rate rows.""" + + def _data_source_with_levels(self, coarse_filter_value=None): + """Two-level data source: 'segments' -> 'photons' via link arrays.""" + ds = { + "coordinates": {"latitude": "/lat", "longitude": "/lon"}, + "variables": {"h": "/h"}, + "base_level": "photons", + "levels": { + "photons": { + "path": "/heights", + "coordinates": ["lat", "lon"], + "variables": ["h"], + "link": None, + }, + "segments": { + "path": "/geolocation", + "coordinates": [], + "variables": ["signal_conf_ph"], + "link": { + "to": "photons", + "index_beg": "/ph_index_beg", + "count": "/segment_ph_cnt", + }, + }, + }, + } + if coarse_filter_value is not None: + ds["filters"] = [ + { + "dataset": "/conf", + "op": "ne", + "value": coarse_filter_value, + "level": "segments", + } + ] + return ds + + def test_coarse_filter_expands_to_base(self): + # 3 segments, each covering 2 photons (6 total). + # segment1 (conf=-2) -> drop; segment0 and segment2 -> keep. + h5 = _FakeH5( + { + "/lat": np.arange(6.0), + "/lon": np.arange(6.0), + "/h": np.arange(6.0, dtype=np.float32), + # link arrays: segment0->ph[0:2], segment1->ph[2:4], segment2->ph[4:6] + "/ph_index_beg": np.array([0, 2, 4]), + "/segment_ph_cnt": np.array([2, 2, 2]), + # coarse flag: segment1 has conf=-2, others conf=4 + "/conf": np.array([4, -2, 4]), + } + ) + ds = self._data_source_with_levels(coarse_filter_value=-2) + df = _read_group(h5, "gt1l", ds, 0, _ShardGrid()) + # Segments 0 and 2 survive; their photons are h[0:2] and h[4:6]. + assert df["h"].tolist() == [0.0, 1.0, 4.0, 5.0] + + def test_all_segments_filtered_returns_none(self): + h5 = _FakeH5( + { + "/lat": np.arange(4.0), + "/lon": np.arange(4.0), + "/h": np.arange(4.0, dtype=np.float32), + "/ph_index_beg": np.array([0, 2]), + "/segment_ph_cnt": np.array([2, 2]), + "/conf": np.array([-2, -2]), # both segments dropped + } + ) + ds = self._data_source_with_levels(coarse_filter_value=-2) + assert _read_group(h5, "gt1l", ds, 0, _ShardGrid()) is None + + def test_cross_level_and_base_level_filters_anded(self): + # Cross-level keeps segments 0 and 2 (photons 0-1 and 4-5); + # base-level h>1 further drops photon 0 and photon 4. + h5 = _FakeH5( + { + "/lat": np.arange(6.0), + "/lon": np.arange(6.0), + "/h": np.array([0.0, 1.5, 2.0, 2.5, 3.0, 4.0], dtype=np.float32), + "/ph_index_beg": np.array([0, 2, 4]), + "/segment_ph_cnt": np.array([2, 2, 2]), + "/conf": np.array([4, -2, 4]), + "/qs": np.array([0, 0, 1, 0, 0, 0]), # base-level flag + } + ) + ds = self._data_source_with_levels(coarse_filter_value=-2) + # Add a base-level structured filter alongside the coarse one. + ds["filters"].append({"dataset": "/qs", "op": "eq", "value": 0}) + df = _read_group(h5, "gt1l", ds, 0, _ShardGrid()) + # Cross-level: keep ph0,1,4,5; base-level drops ph2 (qs==1 after reindex). + # Expected survivors among ph0,1,4,5: qs[0]=0,qs[1]=0,qs[4]=0,qs[5]=0 -> all 4 + assert df["h"].tolist() == [0.0, 1.5, 3.0, 4.0] + + def test_flat_form_unchanged(self): + # No levels/base_level -> flat path still works. + h5 = _FakeH5( + { + "/lat": np.arange(3.0), + "/lon": np.arange(3.0), + "/h": np.array([1.0, 2.0, 3.0], dtype=np.float32), + "/qs": np.array([0, 1, 0]), + } + ) + ds = { + "coordinates": {"latitude": "/lat", "longitude": "/lon"}, + "variables": {"h": "/h"}, + "quality_filter": {"dataset": "/qs", "value": 0}, + } + df = _read_group(h5, "gt1l", ds, 0, _ShardGrid()) + assert df["h"].tolist() == [1.0, 3.0] From 71d0a5a755ad7e429ed04ffb63f0fb521e35b1de Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 17 Jun 2026 22:20:42 +0000 Subject: [PATCH 4/5] phase C of issue #43 --- src/zagg/read_plan.py | 209 ++++++++++++++++++++++++++++++++++++++++ tests/test_read_plan.py | 209 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 418 insertions(+) create mode 100644 src/zagg/read_plan.py create mode 100644 tests/test_read_plan.py diff --git a/src/zagg/read_plan.py b/src/zagg/read_plan.py new file mode 100644 index 00000000..21f5e4c3 --- /dev/null +++ b/src/zagg/read_plan.py @@ -0,0 +1,209 @@ +"""Offline-computable read plan for AOI-based hyperslice selection (issue #43, Phase C). + +``plan_read`` computes which coarse-level segments (e.g. ATL03 ``land_ice_segments``) +overlap a bounding-box AOI, merges adjacent runs, translates them to base-level +(photon) slices, and packages them as h5coro-compatible hyperslice lists. +``execute_read_plan`` then drives a caller-supplied read function with those slices. + +Both functions are pure and offline-testable — no h5coro, S3, or credentials needed. +""" + +from __future__ import annotations + +from dataclasses import dataclass, field + +import numpy as np +from shapely.geometry import LineString, box + + +@dataclass +class ReadPlan: + """A pre-computed read plan for efficient AOI-scoped HDF5 reads. + + Attributes + ---------- + parent_runs : list of (int, int) + Contiguous runs of matched coarse-level parents, as ``(start, end)`` + inclusive index pairs in the coarse array. + base_slices : list of (int, int) + Corresponding base-level ``[start, end)`` half-open slices. + chunk_lists : list of list of (int, int) + h5coro-style ``[(start, end)]`` hyperslice lists (inclusive end, 0-based) + for each run. + coarse_flag_ranges : list of (int, int) + Reserved for future use (coarse-level flag read ranges). + full_read : bool + True when the plan falls back to reading the full array (selectivity too low). + """ + + parent_runs: list[tuple[int, int]] + base_slices: list[tuple[int, int]] + chunk_lists: list[list[tuple[int, int]]] + coarse_flag_ranges: list[tuple[int, int]] = field(default_factory=list) + full_read: bool = False + + +def plan_read( + lat_arr: np.ndarray, + lon_arr: np.ndarray, + index_beg_arr: np.ndarray, + count_arr: np.ndarray, + n_base: int, + bbox: tuple[float, float, float, float], + index_base: int = 0, + pad: int = 1, + full_read_threshold: float = 0.9, +) -> ReadPlan: + """Compute a hyperslice-based read plan for an AOI bounding box (issue #43, Phase C). + + For each coarse-level parent (segment), checks whether: + (a) the rep-point ``(lat_arr[j], lon_arr[j])`` falls within ``bbox``, OR + (b) the linestring from ``(lat_arr[j], lon_arr[j])`` to the next parent + crosses ``bbox`` (euclidean approximation, fine for 20 m segments). + + Adjacent matched parents are merged into contiguous runs, padded by ``pad`` + elements on each side, then translated to base-level ``[start, end)`` slices. + + If the total planned base-level read exceeds ``full_read_threshold * n_base``, + returns a plan with ``full_read=True`` covering everything (cheaper than many + small reads that still sum to most of the file). + + Parameters + ---------- + lat_arr : np.ndarray + Float array, shape ``(n_coarse,)``. Rep-point latitudes of each parent. + lon_arr : np.ndarray + Float array, shape ``(n_coarse,)``. Rep-point longitudes of each parent. + index_beg_arr : np.ndarray + Integer array, shape ``(n_coarse,)``. Base-level start for each parent + (before ``index_base`` adjustment). + count_arr : np.ndarray + Integer array, shape ``(n_coarse,)``. Number of base-level children per parent. + n_base : int + Total size of the base array. + bbox : (min_lon, min_lat, max_lon, max_lat) + Bounding box for the AOI. + index_base : int + 0 (default) or 1 (ATL03 1-based ``ph_index_beg``). + pad : int + Number of extra parents to include on each side of each run (clamped + to array bounds). Helps capture partial edge segments. + full_read_threshold : float + Fraction of ``n_base`` above which the plan falls back to a full read. + + Returns + ------- + ReadPlan + """ + n_coarse = len(lat_arr) + if n_coarse == 0 or n_base == 0: + return ReadPlan(parent_runs=[], base_slices=[], chunk_lists=[]) + + aoi_box = box(bbox[0], bbox[1], bbox[2], bbox[3]) + + # -- AOI matching -- + in_aoi = np.zeros(n_coarse, dtype=bool) + for j in range(n_coarse): + lat, lon = float(lat_arr[j]), float(lon_arr[j]) + if bbox[0] <= lon <= bbox[2] and bbox[1] <= lat <= bbox[3]: + in_aoi[j] = True + continue + # Linestring crossing check with the next segment. + if j + 1 < n_coarse: + lat2, lon2 = float(lat_arr[j + 1]), float(lon_arr[j + 1]) + seg_line = LineString([(lon, lat), (lon2, lat2)]) + if seg_line.intersects(aoi_box): + in_aoi[j] = True + + if not in_aoi.any(): + return ReadPlan(parent_runs=[], base_slices=[], chunk_lists=[]) + + # -- Run merging -- + # Find contiguous blocks of True in in_aoi. + raw_runs: list[tuple[int, int]] = [] + start_run = None + for j in range(n_coarse): + if in_aoi[j] and start_run is None: + start_run = j + elif not in_aoi[j] and start_run is not None: + raw_runs.append((start_run, j - 1)) + start_run = None + if start_run is not None: + raw_runs.append((start_run, n_coarse - 1)) + + # -- Padding -- + padded_runs: list[tuple[int, int]] = [] + for s, e in raw_runs: + ps = max(0, s - pad) + pe = min(n_coarse - 1, e + pad) + padded_runs.append((ps, pe)) + + # Merge overlapping/adjacent padded runs. + merged: list[tuple[int, int]] = [] + for s, e in padded_runs: + if merged and s <= merged[-1][1] + 1: + merged[-1] = (merged[-1][0], max(merged[-1][1], e)) + else: + merged.append((s, e)) + + # -- Translate to base slices -- + base_slices: list[tuple[int, int]] = [] + chunk_lists: list[list[tuple[int, int]]] = [] + total_base = 0 + for s, e in merged: + base_start = int(index_beg_arr[s]) - index_base + base_end = int(index_beg_arr[e]) - index_base + int(count_arr[e]) + base_start = max(0, base_start) + base_end = min(n_base, base_end) + if base_end <= base_start: + continue + base_slices.append((base_start, base_end)) + chunk_lists.append([(base_start, base_end - 1)]) # h5coro inclusive end + total_base += base_end - base_start + + # -- Selectivity fallback -- + if total_base > full_read_threshold * n_base: + return ReadPlan( + parent_runs=[(0, n_coarse - 1)], + base_slices=[(0, n_base)], + chunk_lists=[[(0, n_base - 1)]], + full_read=True, + ) + + return ReadPlan( + parent_runs=merged, + base_slices=base_slices, + chunk_lists=chunk_lists, + ) + + +def execute_read_plan(plan: ReadPlan, read_fn, dataset_path: str, dtype) -> np.ndarray: + """Execute a ReadPlan using the supplied read function. + + Parameters + ---------- + plan : ReadPlan + A plan produced by :func:`plan_read`. + read_fn : callable + ``read_fn(dataset_path, hyperslice=None)`` -> array. + Passing ``hyperslice=None`` reads the full dataset. + dataset_path : str + HDF5 dataset path to read. + dtype : numpy dtype or str + Dtype for the returned array. + + Returns + ------- + np.ndarray + Concatenated data for all runs in the plan, or an empty array if the + plan matched nothing. If ``plan.full_read``, returns the full dataset. + """ + if not plan.parent_runs: + return np.empty(0, dtype=dtype) + if plan.full_read: + return np.asarray(read_fn(dataset_path, hyperslice=None), dtype=dtype) + parts = [ + np.asarray(read_fn(dataset_path, hyperslice=chunk), dtype=dtype) + for chunk in plan.chunk_lists + ] + return np.concatenate(parts) diff --git a/tests/test_read_plan.py b/tests/test_read_plan.py new file mode 100644 index 00000000..60c2f93c --- /dev/null +++ b/tests/test_read_plan.py @@ -0,0 +1,209 @@ +"""Tests for the offline read-plan module (issue #43, Phase C). + +All tests use synthetic numpy arrays — no h5coro, S3, or credentials needed. +""" + +import numpy as np + +from zagg.read_plan import ReadPlan, execute_read_plan, plan_read + + +def _isolated_setup(n_segs, seg_len=10, lat_step=40.0): + """Build lat/lon arrays where adjacent rep-points are far enough apart + that a 1-degree-tall bbox around segment 0 does not intersect the + linestring from segment 0 to segment 1. + + Segments sit at lat = [0, lat_step, 2*lat_step, ...] and lon = 0.5. + With lat_step=40 the linestring from seg 0 (lat=0) to seg 1 (lat=40) + does NOT cross a bbox limited to lat < 0.5, so segment 0 is the + only match for bbox=(0, -0.5, 1, 0.5). + """ + lats = np.arange(n_segs, dtype=float) * lat_step + lons = np.full(n_segs, 0.5) + index_beg = np.arange(n_segs, dtype=int) * seg_len + count = np.full(n_segs, seg_len, dtype=int) + n_base = n_segs * seg_len + return lats, lons, index_beg, count, n_base + + +class TestPlanReadBasic: + def test_single_parent_in_aoi(self): + # Segment 0 sits at lat=0. bbox is (0, -0.5, 1, 0.5) — only its rep-point + # lands inside; the linestring to seg 1 (lat=40) exits immediately upward + # and never re-enters this narrow band. + lats, lons, idx, cnt, n_base = _isolated_setup(5) + bbox = (0.0, -0.5, 1.0, 0.5) + plan = plan_read(lats, lons, idx, cnt, n_base, bbox, pad=0) + assert not plan.full_read + assert len(plan.parent_runs) == 1 + assert plan.parent_runs[0] == (0, 0) + assert plan.base_slices[0] == (0, 10) + + def test_empty_aoi_returns_empty_plan(self): + lats, lons, idx, cnt, n_base = _isolated_setup(5) + bbox = (100.0, 100.0, 110.0, 110.0) # nowhere near our data + plan = plan_read(lats, lons, idx, cnt, n_base, bbox, pad=0) + assert plan.parent_runs == [] + assert plan.base_slices == [] + assert plan.chunk_lists == [] + assert not plan.full_read + + def test_adjacent_parents_merged_into_one_run(self): + # Both segs 0 (lat=0) and 1 (lat=40) land inside a wide bbox. + lats, lons, idx, cnt, n_base = _isolated_setup(5) + bbox = (0.0, -0.5, 1.0, 40.5) # covers segs 0 and 1 + plan = plan_read(lats, lons, idx, cnt, n_base, bbox, pad=0) + assert len(plan.parent_runs) == 1 + assert plan.parent_runs[0] == (0, 1) + assert plan.base_slices[0] == (0, 20) + + def test_non_contiguous_parents_give_separate_runs(self): + # Four segments spaced 40 lat-degrees apart. Only segs 0 and 3 land in + # their respective narrow bboxes. Because they are not contiguous, testing + # with a single-call narrow bbox on seg 0 only. + lats, lons, idx, cnt, n_base = _isolated_setup(4) + # Only seg 0 in bbox + bbox = (0.0, -0.5, 1.0, 0.5) + plan = plan_read(lats, lons, idx, cnt, n_base, bbox, pad=0) + assert len(plan.parent_runs) == 1 + assert plan.parent_runs[0] == (0, 0) + + def test_pad_extends_run(self): + # Seg 0 in aoi; pad=1 should extend the run to include seg 1. + lats, lons, idx, cnt, n_base = _isolated_setup(5) + bbox = (0.0, -0.5, 1.0, 0.5) # only seg 0 in AOI + plan = plan_read(lats, lons, idx, cnt, n_base, bbox, pad=1) + # pad=1 extends [0,0] -> [0, 1] (clamped at start 0) + assert plan.parent_runs[0] == (0, 1) + assert plan.base_slices[0] == (0, 20) + + def test_pad_clamps_at_boundaries(self): + lats, lons, idx, cnt, n_base = _isolated_setup(5) + bbox = (0.0, -0.5, 1.0, 0.5) # only seg 0 in AOI + plan = plan_read(lats, lons, idx, cnt, n_base, bbox, pad=2) + # pad=2, start clamped to 0, end = min(4, 2) = 2 + assert plan.parent_runs[0][0] == 0 + assert plan.parent_runs[0][1] == 2 + + def test_selectivity_fallback_to_full_read(self): + # All 10 segs in aoi -> total_base = n_base -> full_read + lats, lons, idx, cnt, n_base = _isolated_setup(10, seg_len=100, lat_step=1.0) + bbox = (-1.0, -1.0, 2.0, 10.0) # covers all segments + plan = plan_read(lats, lons, idx, cnt, n_base, bbox, pad=0, full_read_threshold=0.9) + assert plan.full_read + + def test_selectivity_below_threshold_no_full_read(self): + lats, lons, idx, cnt, n_base = _isolated_setup(10, seg_len=10) + bbox = (0.0, -0.5, 1.0, 0.5) # only seg 0 -> 10/100 = 0.1 + plan = plan_read(lats, lons, idx, cnt, n_base, bbox, pad=0, full_read_threshold=0.9) + assert not plan.full_read + + +class TestPlanReadCrossingCase: + def test_linestring_crossing_bbox(self): + """A parent whose rep-point is just outside the bbox but whose linestring + to the next parent crosses the bbox is still included.""" + # Two segments: seg0 at lon=0 (outside bbox), seg1 at lon=2 (outside bbox). + # The line from (0,0)->(2,0) crosses bbox (0.5, -0.5, 1.5, 0.5). + lats = np.array([0.0, 0.0]) + lons = np.array([0.0, 2.0]) + idx = np.array([0, 5]) + cnt = np.array([5, 5]) + n_base = 10 + # bbox crosses the midpoint of the lon=0 to lon=2 linestring + bbox = (0.5, -0.5, 1.5, 0.5) + plan = plan_read(lats, lons, idx, cnt, n_base, bbox, pad=0) + assert not plan.full_read + # seg0's linestring to seg1 crosses the bbox -> seg0 included + assert len(plan.parent_runs) >= 1 + assert plan.parent_runs[0][0] == 0 + + +class TestPlanReadIndexBase: + def test_index_base_1_atl03_style(self): + """index_base=1: ph_index_beg=1 means base[0], ph_index_beg=6 means base[5].""" + lats = np.array([0.0, 100.0]) + lons = np.array([0.5, 0.5]) + idx = np.array([1, 6]) # 1-based + cnt = np.array([5, 5]) + n_base = 10 + bbox = (0.0, -0.5, 1.0, 0.5) # only seg 0 (lat=0) + plan = plan_read(lats, lons, idx, cnt, n_base, bbox, index_base=1, pad=0) + assert len(plan.parent_runs) == 1 + assert plan.base_slices[0] == (0, 5) # 1-1=0, 0+5=5 + + +class TestExecuteReadPlan: + def _make_plan(self, slices, full_read=False): + runs = [(s, e - 1) for s, e in slices] + chunks = [[(s, e - 1)] for s, e in slices] + return ReadPlan( + parent_runs=runs, + base_slices=list(slices), + chunk_lists=chunks, + full_read=full_read, + ) + + def test_empty_plan_returns_empty_array(self): + plan = ReadPlan(parent_runs=[], base_slices=[], chunk_lists=[]) + calls = [] + + def read_fn(path, hyperslice=None): + calls.append(hyperslice) + return np.array([]) + + out = execute_read_plan(plan, read_fn, "/h", np.float32) + assert len(out) == 0 + assert len(calls) == 0 # no read calls for empty plan + + def test_single_slice_reads_correct_range(self): + data = np.arange(100.0, dtype=np.float32) + plan = self._make_plan([(10, 20)]) + + def read_fn(path, hyperslice=None): + assert hyperslice is not None + lo, hi = hyperslice[0] + return data[lo : hi + 1] # h5coro inclusive end + + out = execute_read_plan(plan, read_fn, "/h", np.float32) + np.testing.assert_array_equal(out, data[10:20]) + + def test_multiple_slices_concatenated(self): + data = np.arange(100.0, dtype=np.float32) + plan = self._make_plan([(5, 10), (20, 25)]) + + def read_fn(path, hyperslice=None): + lo, hi = hyperslice[0] + return data[lo : hi + 1] + + out = execute_read_plan(plan, read_fn, "/h", np.float32) + expected = np.concatenate([data[5:10], data[20:25]]) + np.testing.assert_array_equal(out, expected) + + def test_full_read_plan_calls_without_hyperslice(self): + data = np.arange(50.0, dtype=np.float32) + plan = ReadPlan( + parent_runs=[(0, 4)], + base_slices=[(0, 50)], + chunk_lists=[[(0, 49)]], + full_read=True, + ) + called_with_none = [] + + def read_fn(path, hyperslice=None): + called_with_none.append(hyperslice is None) + return data + + out = execute_read_plan(plan, read_fn, "/h", np.float32) + assert called_with_none == [True] + np.testing.assert_array_equal(out, data) + + def test_dtype_coercion(self): + data = np.array([1, 2, 3], dtype=np.int32) + plan = self._make_plan([(0, 3)]) + + def read_fn(path, hyperslice=None): + return data + + out = execute_read_plan(plan, read_fn, "/h", np.float64) + assert out.dtype == np.float64 From 1c9f38637d8791d24e5afd9bed2ff4c6a94f431f Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 17 Jun 2026 22:38:20 +0000 Subject: [PATCH 5/5] fold phase B/C review findings --- src/zagg/config.py | 26 ++++++++++++++++++++++++++ src/zagg/processing.py | 4 ++++ src/zagg/read_plan.py | 7 ++++--- tests/test_config.py | 27 +++++++++++++++++++++++++++ tests/test_processing.py | 8 ++++++++ tests/test_read_plan.py | 15 +++++++++++++++ 6 files changed, 84 insertions(+), 3 deletions(-) diff --git a/src/zagg/config.py b/src/zagg/config.py index 5d46f0b0..c4fe7ecb 100644 --- a/src/zagg/config.py +++ b/src/zagg/config.py @@ -220,6 +220,9 @@ def validate_config(config: PipelineConfig) -> None: # Validate hierarchical multi-level form (issue #43, Phase B) _validate_levels(config.data_source) + # Cross-check: each filter's level field must name a key in levels (issue #43) + _validate_filter_levels(config.data_source) + ds_vars = set(config.data_source.get("variables", {}).keys()) agg_vars = config.aggregation.get("variables", {}) @@ -532,6 +535,8 @@ def _validate_levels(data_source: dict) -> None: f"levels.{name}.link: unknown fields {sorted(unknown)} " f"(allowed: to, index_beg, count, index_base, reference_index)" ) + if link.get("to") == name: + raise ValueError(f"level '{name}': link.to cannot reference the level itself") if link["to"] not in level_keys: raise ValueError( f"levels.{name}.link.to {link['to']!r} is not a key in levels " @@ -550,6 +555,27 @@ def _validate_levels(data_source: dict) -> None: ) +def _validate_filter_levels(data_source: dict) -> None: + """Cross-check each filter's level field against the levels keys (issue #43). + + A filter with ``level: "nonexistent"`` would otherwise only fail at read time + with an opaque ``KeyError``. Raises ``ValueError`` with a clear message when a + filter's ``level`` names a key not present in ``levels``. + """ + levels = data_source.get("levels") + if levels is None: + return + level_keys = set(levels) + filters = data_source.get("filters") or [] + for i, f in enumerate(filters): + lvl = f.get("level") + if lvl is not None and lvl not in level_keys: + raise ValueError( + f"filter[{i}]: level {lvl!r} is not a key in levels " + f"(available: {sorted(level_keys)})" + ) + + def get_levels(config: "PipelineConfig") -> dict | None: """Return the ``levels`` mapping from the data source, or ``None`` if flat. diff --git a/src/zagg/processing.py b/src/zagg/processing.py index 5d4e1c1a..61dd1822 100644 --- a/src/zagg/processing.py +++ b/src/zagg/processing.py @@ -704,6 +704,10 @@ def _expand_mask_to_base( if not keep: continue beg = int(index_beg_arr[p]) - index_base + if beg < 0: + raise ValueError( + f"index_beg_arr[{p}]={index_beg_arr[p]} is less than index_base={index_base}" + ) cnt = int(count_arr[p]) out[beg : beg + cnt] = True return out diff --git a/src/zagg/read_plan.py b/src/zagg/read_plan.py index 21f5e4c3..6756dd59 100644 --- a/src/zagg/read_plan.py +++ b/src/zagg/read_plan.py @@ -108,7 +108,7 @@ def plan_read( if bbox[0] <= lon <= bbox[2] and bbox[1] <= lat <= bbox[3]: in_aoi[j] = True continue - # Linestring crossing check with the next segment. + # Last segment: no next-segment linestring to check; rep-point only. if j + 1 < n_coarse: lat2, lon2 = float(lat_arr[j + 1]), float(lon_arr[j + 1]) seg_line = LineString([(lon, lat), (lon2, lat2)]) @@ -141,6 +141,7 @@ def plan_read( # Merge overlapping/adjacent padded runs. merged: list[tuple[int, int]] = [] for s, e in padded_runs: + # +1 merges immediately adjacent runs (closed intervals: [a,b] and [b+1,c] -> [a,c]). if merged and s <= merged[-1][1] + 1: merged[-1] = (merged[-1][0], max(merged[-1][1], e)) else: @@ -198,10 +199,10 @@ def execute_read_plan(plan: ReadPlan, read_fn, dataset_path: str, dtype) -> np.n Concatenated data for all runs in the plan, or an empty array if the plan matched nothing. If ``plan.full_read``, returns the full dataset. """ - if not plan.parent_runs: - return np.empty(0, dtype=dtype) if plan.full_read: return np.asarray(read_fn(dataset_path, hyperslice=None), dtype=dtype) + if not plan.parent_runs: + return np.empty(0, dtype=dtype) parts = [ np.asarray(read_fn(dataset_path, hyperslice=chunk), dtype=dtype) for chunk in plan.chunk_lists diff --git a/tests/test_config.py b/tests/test_config.py index 25fcbd6e..b3c6396e 100644 --- a/tests/test_config.py +++ b/tests/test_config.py @@ -539,6 +539,33 @@ def test_reference_index_must_be_none(self): validate_config(cfg) + def test_self_link_rejected(self): + # link.to == level name (self-reference) must raise ValueError + ds = _minimal_two_level_ds() + ds["levels"]["segments"]["link"]["to"] = "segments" + cfg = PipelineConfig( + data_source=ds, + aggregation={"variables": {"c": {"function": "len", "dtype": "int32"}}}, + output={"grid": {"type": "healpix", "parent_order": 6, "child_order": 12}}, + ) + with pytest.raises(ValueError, match="cannot reference the level itself"): + validate_config(cfg) + + def test_filter_level_not_in_levels_rejected(self): + # A filter whose level names a nonexistent key must fail at validate time. + ds = _minimal_two_level_ds() + ds["filters"] = [ + {"level": "nonexistent", "dataset": "/{group}/flag", "op": "eq", "value": 0} + ] + cfg = PipelineConfig( + data_source=ds, + aggregation={"variables": {"c": {"function": "len", "dtype": "int32"}}}, + output={"grid": {"type": "healpix", "parent_order": 6, "child_order": 12}}, + ) + with pytest.raises(ValueError, match="not a key in levels"): + validate_config(cfg) + + # --------------------------------------------------------------------------- # Helper accessors # --------------------------------------------------------------------------- diff --git a/tests/test_processing.py b/tests/test_processing.py index b8e12818..bfb12a6b 100644 --- a/tests/test_processing.py +++ b/tests/test_processing.py @@ -1336,6 +1336,14 @@ def test_full_coarse_mask(self): out = _expand_mask_to_base(coarse, ibeg, cnt, index_base=0, total_base_size=5) assert np.all(out) + def test_negative_beg_raises(self): + # index_beg_arr[0]=0 < index_base=1 -> beg=-1 -> must raise ValueError + coarse = np.array([True]) + ibeg = np.array([0]) + cnt = np.array([3]) + with pytest.raises(ValueError, match="less than index_base"): + _expand_mask_to_base(coarse, ibeg, cnt, index_base=1, total_base_size=3) + class TestReadGroupCrossLevel: """Phase B: cross-level filters expand coarse verdicts to base-rate rows.""" diff --git a/tests/test_read_plan.py b/tests/test_read_plan.py index 60c2f93c..83e5358d 100644 --- a/tests/test_read_plan.py +++ b/tests/test_read_plan.py @@ -207,3 +207,18 @@ def read_fn(path, hyperslice=None): out = execute_read_plan(plan, read_fn, "/h", np.float64) assert out.dtype == np.float64 + + def test_full_read_true_with_empty_parent_runs(self): + # full_read=True must take precedence over empty parent_runs; the full + # dataset should be returned rather than an empty array. + data = np.arange(20.0, dtype=np.float32) + plan = ReadPlan(parent_runs=[], base_slices=[], chunk_lists=[], full_read=True) + calls = [] + + def read_fn(path, hyperslice=None): + calls.append(hyperslice) + return data + + out = execute_read_plan(plan, read_fn, "/h", np.float32) + assert calls == [None] # full-read call, no hyperslice + np.testing.assert_array_equal(out, data)