diff --git a/deployment/aws/lambda_handler.py b/deployment/aws/lambda_handler.py index 7857f0e..be04811 100644 --- a/deployment/aws/lambda_handler.py +++ b/deployment/aws/lambda_handler.py @@ -134,13 +134,11 @@ # Import cloud-agnostic processing from zagg.config import get_handoff, get_store_layout, load_config_from_dict -from zagg.grids.base import shard_label from zagg.processing import ( write_dataframe_to_zarr, write_ragged_to_zarr, write_shard_to_zarr, ) -from zagg.processing.write import _block_index_key from zagg.store import open_store # Set up structured logging @@ -935,11 +933,6 @@ def _handle_process(event: Dict[str, Any], context: Any) -> Dict[str, Any]: sharded = getattr(grid, "sharded", False) store_path = event["store_path"] shard_key = event["shard_key"] - # K==1 vs K>1 is fixed by the grid, not the chunk count (issue #91), so the - # streaming callback can pick the ragged key without a materialized list: at - # K==1 the chunk IS the shard (keyed by ``shard_key``); at K>1 each finer - # chunk is keyed by its own block index. - single_chunk = int(getattr(grid, "chunks_per_shard", 1)) == 1 store_box: dict = {} write_error: dict = {} @@ -1017,15 +1010,9 @@ def _write_chunk(block_index, carrier, ragged): # write_dataframe_to_zarr no-ops on an empty carrier, so no per-chunk # emptiness check is needed. Use each chunk's own block_index. write_dataframe_to_zarr(carrier, store, grid=grid, chunk_idx=block_index) - # At K==1 the CSR subgroup is named by the grid's shard label (the - # decimal morton string for HEALPix — issue #199); K>1 keeps the - # flattened block-index int. Mirrors runner._process_and_write. - ragged_key = ( - shard_label(grid, shard_key) - if single_chunk - else _block_index_key(block_index, grid) - ) - write_ragged_to_zarr(ragged, store, grid=grid, shard_key=ragged_key) + # Ragged fields land in their vlen-bytes arrays at the same + # block (issue #209). Mirrors runner._process_and_write. + write_ragged_to_zarr(ragged, store, grid=grid, chunk_idx=block_index) except Exception as e: # Mirror the buffered path's ``except``: record the failure, stop # writing, and let the run surface a 500 after process_shard returns. diff --git a/docs/ragged_layout.md b/docs/ragged_layout.md new file mode 100644 index 0000000..1348d6b --- /dev/null +++ b/docs/ragged_layout.md @@ -0,0 +1,179 @@ +# Ragged store layout (zagg-ragged/1) + +A `kind: ragged` field (a per-cell t-digest, a per-cell photon list — anything +whose per-cell length varies) is stored as **one self-describing +`variable_length_bytes` array** on the cell grid. This is the +[issue #209](https://github.com/englacial/zagg/issues/209) layout; it replaced +the per-inner-chunk CSR subgroups (`values`/`offsets`/`cell_ids`, ~7 objects per +populated inner chunk) with a single vlen array per field. The convention is +versioned `zagg-ragged/1` and stamped into the array's attrs. + +The design goals it meets: **one object per shard** on the write side (deleting +the ~K×7 tiny-PUT storm and the [issue #142](https://github.com/englacial/zagg/issues/142) +write-fanout thread pool that existed only to parallelize it), **2-GET random +access** to any single cell on the read side, and a wire format frozen tightly +enough that the [issue #210](https://github.com/englacial/zagg/issues/210) +typed-dtype migration is metadata-only. + +## Layout + +A ragged field `{field}` under a product group is up to three sibling arrays: + +``` +{group}/{field} <- vlen array; populated cell i holds the raw + little-endian bytes of its (n, *inner_shape) payload +{group}/{field}_locations <- LOCATED fields only (issue #87): the per-row uint64 + location words, row-aligned with {field} +{group}/morton <- per-cell uint64 morton coordinate (zagg's standard + HEALPix coordinate array; the chunk-identity source) +``` + +Each populated cell's value is the raw little-endian bytes of an +`(n, *inner_shape)` array (`n` varies per cell — e.g. `(k_centroids, 2)` for a +t-digest whose `inner_shape` is `(2,)`). Empty cells keep the `b""` fill, and an +inner chunk with no ragged data is omitted from disk entirely — the same +sub-shard sparsity the dense arrays get. + +The array's codec chain is `[vlen-bytes, zstd(level=3)]`. The zstd deviates from +the dense arrays' bytes-only/uncompressed policy deliberately: a vlen payload has +no fixed-width raw layout to preserve, and level 3 matches the coverage-sidecar +precedent (`zagg.hive._ZSTD_LEVEL`), fixed so identical payloads produce +identical objects across workers. + +## The attrs contract + +The element interpretation is **self-describing** in the array's attrs under the +`ragged` key (`grids.base.RAGGED_ELEMENT_ATTR`), so a reader decodes exactly what +the writer declared rather than hardcoding a dtype: + +```json +"ragged": { + "spec": "zagg-ragged/1", + "element": {"dtype": "float32", "shape": [-1, 2]}, + "locations": "h_tdigest_locations" +} +``` + +- **`spec`** (`grids.base.RAGGED_SPEC`) — the convention version. Readers + strict-check it (`readers/tdigest_tensor._open_ragged`): an unknown/future spec + **raises**, never half-parses. This is the coverage-envelope discipline applied + to the ragged layout. +- **`element`** — `{"dtype": "", "shape": [-1, *inner_shape]}`. The + `-1` marks the per-cell varying count, so a reader reconstructs cell `i` as + `np.frombuffer(a[i], dtype).reshape(-1, *inner_shape)`. The bytes are always + little-endian, independent of the producing machine. +- **`locations`** — present only on a LOCATED field's payload array (issue #87); + its value is the name of the sibling uint64 array carrying the per-row location + words. A reader binds the location channel **by this declaration**, never by + reconstructing the `{field}_locations` naming convention. An unlocated field + records nothing here. + +A vlen array without a well-formed `element` declaration is **not** a zagg ragged +array — pre-issue-209 CSR stores are a hard break, and the readers raise a +pointed error rather than decode under a guessed layout. + +## Wire framing (golden-pinned) + +Within one inner chunk the `variable_length_bytes` codec frames the cells before +compression as (little-endian throughout): + +``` +u32 cell_count +per cell: u32 payload_length || payload_bytes +``` + +i.e. numcodecs' `VLenBytes`/`VLenArray` framing — a `u32` count of cells, then +for each cell a `u32` byte length followed by that many payload bytes (`0` for an +empty cell). The `payload_bytes` are `np.ascontiguousarray(value).tobytes()` of +the cell's `(n, *inner_shape)` array in the declared dtype. + +This exact byte vector is frozen by a golden test +(`tests/test_processing.py::TestRaggedVlenLayout::test_golden_inner_chunk_framing`): +round-trip tests pass under any self-consistent encoding, so only a fixed byte +vector pins the convention. It is what guarantees the later metadata-only +migration to a typed `vlen-array` dtype (byte-compatible with numcodecs +`VLenArray`) without rewriting data — see the +[migration note](#issue-210-typed-dtype-migration). + +## Sharded vs per-inner-chunk layouts + +Both layouts hold the same logical data and are self-describing in the array's +own metadata (its `chunk_grid` and whether it has a `shards` outer chunk), so a +single reader code path reads either. Which one a product gets depends on the +write path (`grids.base.ragged_array_spec`, `shard_shape` argument): + +| layout | when | on disk | single-cell read | +|---|---|---|---| +| **sharded** (`ShardingCodec`) | sharded flat path (`write_shard_to_zarr`) and every hive leaf (`write_ragged_leaf_to_zarr`) | ONE object per shard; the shard's K inner chunks live inside it with an internal index | 2 GETs (index suffix + one ranged inner chunk) | +| **per-inner-chunk** (regular array) | UNSHARDED per-chunk write (`write_ragged_to_zarr`, the runner / Lambda streaming callback) | one object per inner chunk | 1 GET (the object) | + +The GET counts are for the data objects only and exclude the one-time array-open +metadata read (amortized across all cells of a store). The sharded 2-GET count is +pinned by `test_two_ranged_gets_on_sharded_store`; the unsharded 1-GET count is +analytic (a regular array indexes the single chunk holding the cell). + +**Why the unsharded flat K>1 path keeps per-inner-chunk objects** (the review's +Q1 resolution): the streaming write path writes each chunk independently as it is +produced, then frees it — the [issue #91](https://github.com/englacial/zagg/issues/91) +stream-and-free bound. Folding all K chunks into one sharded object would force a +read-modify-write of that shared object on every chunk, defeating stream-and-free +and re-introducing the memory the sharded slab pass is careful to bound. So the +regular-chunked (one object per inner chunk) layout is retained there. It is +**not** a reader fork: the reader derives the stored span from `arr.shards or +arr.chunks` and reads either identically (pinned by +`test_sharded_and_regular_layouts_read_identically` — same logical data through +both writers yields identical tensors and chunk ids). + +Empty inner chunks are omitted from the sharded object's index (the `2^64-1` +sentinel in the shard footer of K `(offset, nbytes)` u64 pairs + crc32c), so +object size scales with **populated** chunks only. + +## The hive-leaf reader contract + +A [hive](hive_layout.md) leaf zarr is exactly this layout scoped to one shard. +Under the leaf's `{group}` path a reader finds: + +- the ragged vlen array `{group}/{field}` with its versioned `ragged` attrs, +- the sibling `morton` coordinate array (chunk identity), +- and, for a located field, the `{field}_locations` sibling — + +all sharded as one whole-leaf `ShardingCodec` object (one stored span). So a hive +product is **read one leaf at a time**: open the leaf store +(`hive.shard_leaf_path`) and pass the same `field` path to the readers. The +readers are store-scoped and never traverse the hive digit tree — leaf discovery +is the coverage MOC's / walker's job ([issue #200](https://github.com/englacial/zagg/issues/200)). +The flat-sharded and hive-leaf writers are pinned to store byte-identical +per-cell payloads (`test_hive_leaf_parity_with_flat_sharded`) so the two backends +cannot drift. + +## Read paths + +`zagg.readers.tdigest_tensor` consumes the layout two ways: + +- **Whole-store sweep** (`read_tensors`, `read_raw_values`, `read_locations`) — + one LIST of the array's stored `c/` objects (`_stored_chunk_spans`), + then a per-read-chunk decode. The sweep visits only written data; each stored + object is read in one slice (a sharded object's index suffix is fetched once, + not re-fetched per inner chunk). Each read chunk is a square `(side, side)` + block of cells (`side = isqrt(cells_per_chunk)`, 64 for the production + `chunk_inner` configs), and its coverage-cell morton id is derived from the + sibling `morton` coordinate coarsened to the chunk order. +- **Single-cell random access** (`read_cell`) — indexes the vlen array directly. + On a sharded store that is exactly **2 ranged GETs** (the shard-index suffix, + then the one inner chunk holding the cell), never the whole shard object. An + out-of-range index raises `IndexError` naming the valid range (no negative-index + wrap); an absent cell returns the zero-length `(0, *inner_shape)` array. Works + on the `{field}_locations` sibling too. + +## Issue #210 typed-dtype migration + +The `ragged` attrs block is the **interim** element contract. The +[issue #210](https://github.com/englacial/zagg/issues/210) migration moves the +element declaration (dtype + inner shape) into the zarr data type itself — a +typed `vlen-array` dtype — and supersedes (bumps or removes) the `spec` +marker. Because the on-wire cell framing is already byte-compatible with +numcodecs `VLenArray` and frozen by the golden test above, that migration is +**metadata-only**: the stored chunk bytes do not change, only the array's +declared data type and the `spec` gate. Existing `zagg-ragged/1` data stays +readable; a reader that predates the bump fails loudly on the new `spec` rather +than mis-decoding. diff --git a/mkdocs.yml b/mkdocs.yml index e7a4dde..b14f389 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -26,6 +26,7 @@ nav: - "Architecture": "design/architecture.md" - "Schema": "design/schema.md" - "Sharded output": "sharding.md" + - "Ragged store layout": "ragged_layout.md" - "Hive store layout": "hive_layout.md" - "Typed morton boundary": "morton_arrow.md" - "Deployment": diff --git a/notebooks/tdigest_reader_example.ipynb b/notebooks/tdigest_reader_example.ipynb index 12bfc67..dd89d51 100644 --- a/notebooks/tdigest_reader_example.ipynb +++ b/notebooks/tdigest_reader_example.ipynb @@ -4,7 +4,40 @@ "cell_type": "markdown", "id": "ac52de6c", "metadata": {}, - "source": "# Reading t-digest products: tensors and raw values\n\n[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/englacial/zagg/main?urlpath=lab/tree/notebooks/tdigest_reader_example.ipynb)\n\n_Runs end-to-end on [Binder](https://mybinder.org/v2/gh/englacial/zagg/main?urlpath=lab/tree/notebooks/tdigest_reader_example.ipynb): it builds a small synthetic t-digest store in memory -- no cloud data, no credentials. The Binder image already provides `zagg[analysis]` (incl. matplotlib) via the repo's `.binder/` environment._\n\nzagg writes a gridded Zarr product where each coverage chunk is a `64 x 64`\nblock of child cells, and each populated cell carries a **t-digest** sketch of\nits observed values (e.g. ICESat-2 photon elevations). The sketch is stored as\na ragged [CSR](https://en.wikipedia.org/wiki/Sparse_matrix#Compressed_sparse_row_(CSR,_CRS_or_Yale_format))\nfield (issue #48).\n\nA downstream client usually wants something denser and fixed-size. The\n`zagg.readers` package (issue #79) reconstructs that from the stored digests:\n\n- **`read_tensors`** -- a generator yielding `(tensor, morton_index)` per chunk,\n where `tensor` has shape `(64, 64, n_bins)`: each cell's digest rasterized into\n `n_bins` evenly-spaced z-bins.\n- **`read_raw_values`** -- the lossless companion: when a cell's digest was built\n with no merges, its centroid means *are* the original samples, so the raw value\n vector is recovered exactly (and it raises when that is not possible).\n\nThis notebook is **self-contained**: it builds a small synthetic t-digest store\nin memory, so it runs anywhere (including Binder) with no external data or\ncredentials. The store is built with the same public API zagg uses to write the\nreal product (`zagg.stats.tdigest.build_tdigest` + `zagg.csr.write_csr`), so the\nread path exercised here is identical to the one a real product uses." + "source": [ + "# Reading t-digest products: tensors and raw values\n", + "\n", + "[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/englacial/zagg/main?urlpath=lab/tree/notebooks/tdigest_reader_example.ipynb)\n", + "\n", + "_Runs end-to-end on [Binder](https://mybinder.org/v2/gh/englacial/zagg/main?urlpath=lab/tree/notebooks/tdigest_reader_example.ipynb): it builds a small synthetic t-digest store in memory -- no cloud data, no credentials. The Binder image already provides `zagg[analysis]` (incl. matplotlib) via the repo's `.binder/` environment._\n", + "\n", + "zagg writes a gridded Zarr product where each coverage chunk is a `64 x 64`\n", + "block of child cells, and each populated cell carries a **t-digest** sketch of\n", + "its observed values (e.g. ICESat-2 photon elevations). The sketch is stored in a\n", + "ragged **variable-length-bytes** field (issue #209): one vlen array on the cell\n", + "grid whose populated cells each hold the raw little-endian bytes of their\n", + "`(k_centroids, 2)` `(mean, weight)` digest. The layout is self-describing -- see\n", + "the [ragged layout doc](https://github.com/englacial/zagg/blob/main/docs/ragged_layout.md).\n", + "\n", + "A downstream client usually wants something denser and fixed-size. The\n", + "`zagg.readers` package (issue #79) reconstructs that from the stored digests:\n", + "\n", + "- **`read_tensors`** -- a generator yielding `(tensor, morton_index)` per chunk,\n", + " where `tensor` has shape `(64, 64, n_bins)`: each cell's digest rasterized into\n", + " `n_bins` evenly-spaced z-bins.\n", + "- **`read_raw_values`** -- the lossless companion: when a cell's digest was built\n", + " with no merges, its centroid means *are* the original samples, so the raw value\n", + " vector is recovered exactly (and it raises when that is not possible).\n", + "- **`read_cell`** -- random access to one cell's digest (2 ranged GETs on a\n", + " sharded product, never the whole shard).\n", + "\n", + "This notebook is **self-contained**: it builds a small synthetic t-digest store\n", + "in memory, so it runs anywhere (including Binder) with no external data or\n", + "credentials. The store is built with the same public API zagg uses to write the\n", + "real product (`zagg.stats.tdigest.build_tdigest` +\n", + "`zagg.processing.write_ragged_to_zarr`), so the read path exercised here is\n", + "identical to the one a real product uses." + ] }, { "cell_type": "markdown", @@ -19,11 +52,12 @@ "source": [ "## 1. Build a synthetic t-digest store\n", "\n", - "We mimic two coverage chunks (parent morton ids `100` and `250`), each with a\n", - "handful of populated cells. For every populated cell we draw some samples and\n", - "sketch them with `build_tdigest`, then write all of a chunk's per-cell digests\n", - "as one CSR subgroup under `{field}/{parent_morton}` -- exactly the layout\n", - "`read_tensors` consumes.\n", + "We mimic two coverage chunks at real order-6 morton ids, each with a handful of\n", + "populated cells. For every populated cell we draw some samples and sketch them\n", + "with `build_tdigest`, then write the chunk's per-cell digests into the product's\n", + "ragged vlen array with `write_ragged_to_zarr` -- exactly the layout the readers\n", + "consume. Each chunk's coverage-cell morton id is recovered by the readers from\n", + "the sibling `morton` coordinate array (not from a subgroup name).\n", "\n", "A real product would instead point the readers at an on-disk or S3-backed Zarr\n", "store of the same shape; nothing else about the read code changes." @@ -31,7 +65,7 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "id": "abe4f114", "metadata": { "execution": { @@ -41,41 +75,69 @@ "shell.execute_reply": "2026-06-25T19:05:17.623535Z" } }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "chunk subgroups written: [100, 250]\n" - ] - } - ], + "outputs": [], "source": [ "import numpy as np\n", + "import zarr\n", "from zarr.storage import MemoryStore\n", "\n", - "from zagg.csr import write_csr\n", + "from zagg.config import PipelineConfig\n", + "from zagg.grids import HealpixGrid\n", + "from zagg.grids.morton import morton_word\n", + "from zagg.processing import write_ragged_to_zarr\n", "from zagg.stats.tdigest import build_tdigest\n", "\n", - "FIELD = \"h_tdigest\" # the t-digest field name, e.g. elevation digests\n", + "# One ragged t-digest field on a HEALPix grid: order-6 shards (parent), order-12\n", + "# cells (child) -> each 4096-cell coverage chunk is a 64x64 child block.\n", + "config = PipelineConfig(\n", + " data_source={\"groups\": [\"g\"]},\n", + " aggregation={\n", + " \"coordinates\": {\"morton\": {\"dtype\": \"uint64\", \"fill_value\": 0}},\n", + " \"variables\": {\n", + " \"h_tdigest\": {\n", + " \"function\": \"zagg.stats.tdigest.build_tdigest\",\n", + " \"source\": \"h\",\n", + " \"kind\": \"ragged\",\n", + " \"inner_shape\": [2], # each centroid is a (mean, weight) pair\n", + " \"dtype\": \"float32\",\n", + " \"fill_value\": 0,\n", + " }\n", + " },\n", + " },\n", + " output={\"grid\": {\"type\": \"healpix\", \"parent_order\": 6, \"child_order\": 12}},\n", + ")\n", + "grid = HealpixGrid(6, 12, layout=\"fullsphere\", config=config)\n", + "FIELD = f\"{grid.group_path}/h_tdigest\" # array path the readers consume, e.g. \"12/h_tdigest\"\n", "rng = np.random.default_rng(0)\n", "\n", "\n", - "def write_chunk(store, field, parent_morton, cell_to_values, *, delta=512):\n", - " \"\"\"Write one chunk's per-cell digests as a CSR subgroup {field}/{morton}.\"\"\"\n", + "def write_chunk(store, morton_key, cell_to_values, *, delta=512):\n", + " \"\"\"Sketch each cell's samples and write the chunk into the ragged vlen array.\n", + "\n", + " Returns ``(morton_word, cell_base)`` -- the chunk's packed coverage-cell id\n", + " and the global offset of its first cell (so single cells can be addressed by\n", + " ``read_cell`` below).\n", + " \"\"\"\n", + " word = morton_word(morton_key)\n", + " block = grid.block_index(word)\n", + " base = block[0] * grid.cells_per_chunk\n", + " # The dense per-cell morton coordinate the readers use for chunk identity.\n", + " morton_arr = zarr.open_array(store, path=f\"{grid.group_path}/morton\", mode=\"r+\")\n", + " morton_arr[base : base + grid.cells_per_chunk] = grid.children(word)\n", " cell_ids = sorted(cell_to_values)\n", " digests = [build_tdigest(np.asarray(cell_to_values[c]), delta=delta) for c in cell_ids]\n", - " write_csr(store, f\"{field}/{parent_morton}\", digests, cell_ids)\n", + " write_ragged_to_zarr({\"h_tdigest\": (digests, cell_ids)}, store, grid=grid, chunk_idx=block)\n", + " return word, base\n", "\n", "\n", "store = MemoryStore()\n", + "grid.emit_template(store)\n", "\n", - "# Chunk 100: elevations clustered around ~20 m. cell_id is the cell's position\n", - "# in the chunk's row-major (64x64) children block (so 0 -> (0,0), 4095 -> (63,63)).\n", - "write_chunk(\n", + "# Chunk A: elevations clustered around ~20 m. cell_id is the cell's position in\n", + "# the chunk's row-major (64x64) children block (so 0 -> (0,0), 4095 -> (63,63)).\n", + "word_a, base_a = write_chunk(\n", " store,\n", - " FIELD,\n", - " 100,\n", + " \"1121121\",\n", " {\n", " 0: rng.normal(20.0, 2.0, 3_000),\n", " 5: rng.normal(22.0, 1.5, 2_000),\n", @@ -83,21 +145,17 @@ " },\n", ")\n", "\n", - "# Chunk 250: a different elevation band, ~50 m.\n", - "write_chunk(\n", + "# Chunk B: a different elevation band, ~50 m.\n", + "word_b, base_b = write_chunk(\n", " store,\n", - " FIELD,\n", - " 250,\n", + " \"2431123\",\n", " {\n", " 7: rng.normal(50.0, 3.0, 2_500),\n", " 63: rng.normal(52.0, 2.0, 2_000),\n", " },\n", ")\n", "\n", - "import zarr\n", - "\n", - "group = zarr.open_group(store, path=FIELD, mode=\"r\")\n", - "print(\"chunk subgroups written:\", sorted(int(k) for k in group.group_keys()))\n" + "print(\"chunk morton ids written:\", sorted((word_a, word_b)))" ] }, { @@ -105,14 +163,17 @@ "id": "0c6d98a4", "metadata": {}, "source": [ - "Each chunk is now three CSR arrays per subgroup: `values` (the flat\n", - "concatenation of every cell's `(k_centroids, 2)` `(mean, weight)` digest),\n", - "`offsets`, and `cell_ids`. We can peek at one cell's digest to confirm." + "The field is now ONE `variable_length_bytes` array on the cell grid: each\n", + "populated cell holds the raw little-endian bytes of its `(k_centroids, 2)`\n", + "`(mean, weight)` digest, and empty cells keep the `b\"\"` fill. The element\n", + "interpretation (dtype + shape) rides in the array's `ragged` attrs, so a reader\n", + "decodes exactly what the writer declared. We can peek at one cell's digest with\n", + "`read_cell`, which random-accesses a single cell." ] }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "id": "ee99c084", "metadata": { "execution": { @@ -122,30 +183,13 @@ "shell.execute_reply": "2026-06-25T19:05:17.636536Z" } }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "chunk 100 has 3 populated cells: [0, 5, 4095]\n", - "\n", - "cell 0 digest shape (554, 2) (k_centroids, [mean, weight]):\n", - "[[12.201157 1. ]\n", - " [12.45545 1. ]\n", - " [13.6053095 1. ]\n", - " [13.770982 1. ]\n", - " [13.787327 1. ]]\n" - ] - } - ], + "outputs": [], "source": [ - "from zagg.csr import iter_csr_cells, read_csr\n", + "from zagg.readers import read_cell\n", "\n", - "csr = read_csr(store, f\"{FIELD}/100\")\n", - "cells = iter_csr_cells(csr)\n", - "print(f\"chunk 100 has {len(cells)} populated cells: {[c for c, _ in cells]}\")\n", - "cell_id, digest = cells[0]\n", - "print(f\"\\ncell {cell_id} digest shape {digest.shape} (k_centroids, [mean, weight]):\")\n", + "# read_cell addresses ONE cell by its global position (chunk base + cell id).\n", + "digest = read_cell(store, FIELD, base_a + 5) # cell 5 of chunk A\n", + "print(f\"cell 5 of chunk A: digest shape {digest.shape} (k_centroids, [mean, weight]):\")\n", "print(digest[:5])" ] }, @@ -170,7 +214,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": null, "id": "b7c17ee2", "metadata": { "execution": { @@ -180,16 +224,7 @@ "shell.execute_reply": "2026-06-25T19:05:17.658853Z" } }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "chunk morton=100: tensor (64, 64, 128) dtype=uint32, 3 populated cells, total counts=6404\n", - "chunk morton=250: tensor (64, 64, 128) dtype=uint32, 2 populated cells, total counts=4384\n" - ] - } - ], + "outputs": [], "source": [ "from zagg.readers import read_tensors\n", "\n", @@ -207,16 +242,11 @@ "cell_type": "markdown", "id": "79a77921", "metadata": {}, - "source": [ - "The morton id is the chunk's coverage cell, recovered straight from the store\n", - "(the CSR subgroup name, or a sibling `{field}/morton` coordinate array when\n", - "present). Populated cells land at their row-major position; everything else is\n", - "zero." - ] + "source": "The morton id is the chunk's coverage cell, recovered from the sibling\n`morton` coordinate array (`{group}/morton`, alongside the field under the same\ngroup — the per-cell morton word coarsened to the chunk's order). Populated\ncells land at their row-major position; everything else is zero." }, { "cell_type": "code", - "execution_count": 4, + "execution_count": null, "id": "fa454469", "metadata": { "execution": { @@ -226,19 +256,9 @@ "shell.execute_reply": "2026-06-25T19:05:17.665324Z" } }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "cell 5 (0, 5) count sum: 2000\n", - "cell 4095 (63,63) count sum: 1424\n", - "cell 10 (0,10) count sum: 0 (unpopulated)\n" - ] - } - ], + "outputs": [], "source": [ - "tensor = tensors[100]\n", + "tensor = tensors[word_a]\n", "# cell 5 -> (row 0, col 5); cell 4095 -> (row 63, col 63); cell 10 unpopulated.\n", "print(\"cell 5 (0, 5) count sum:\", int(tensor[0, 5].sum()))\n", "print(\"cell 4095 (63,63) count sum:\", int(tensor[63, 63].sum()))\n", @@ -259,7 +279,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": null, "id": "2029227c", "metadata": { "execution": { @@ -269,29 +289,19 @@ "shell.execute_reply": "2026-06-25T19:05:18.270382Z" } }, - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAArIAAAFUCAYAAADYjN+CAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlcelbwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAP2RJREFUeJzt3Xd0VHX+//HXJCGVFEIJLUR6l4B0QRdCk46FroBfcEEQKSoguOAigh1UVlAEXUVAdDUgQRREICgsLE0B6SWETkghJCHl/v7wZH4Ok4QJTDK54fk4Z85hPvczd9733hnmlc987h2LYRiGAAAAAJNxc3UBAAAAwO0gyAIAAMCUCLIAAAAwJYIsAAAATIkgCwAAAFMiyAIAAMCUCLIAAAAwJYIsAAAATIkgCwAAAFMiyOKWTp48KYvFokWLFuXZ5mrx8fEaMmSIQkJCZLFYNGbMmNte1549e2SxWLR8+XInVmi+GoCC8Mknn8hiseiPP/7Is99zzz0nDw8Ph9b5yiuvyGKxKDU11Rkl3tXeeecdValSRWlpaU5fN8fp1mbNmqVatWopPT3d1aWYAkEWBa58+fKyWCx2t969ezv1eWbMmKG1a9dq69atMgxD77//vlPXbzanTp3SG2+8oWbNmslisahly5Z59l+8eLEaNmwob29vVaxYUWPHjlViYuId90X+TJs2TRaLRRkZGdRRCO6W7XTUlStXNHPmTE2bNk1eXl6uLscpDMPQli1bNGbMGFWoUEEWiyXXz4f169fn+HllsVj01VdfFUq9zz77rK5evXrXf4Y5iiCLQjFo0CAZhmFz+/bbb536HJs2bVKrVq1Uo0YNp67XrP7v//5PFy9e1L/+9S9VqlQpz75z587ViBEjrP+BrlmzRmvWrFGPHj2UlZV1232BwjJt2jQZhiFvb29Xl2JqH3zwgbKysvTEE0+4uhSn2bFjh6ZOnao6deronXfecegxa9eutfvMevTRRwu40j+VLFlSQ4YM0RtvvMEfWA4gyKLYuHz5snx8fFxdRpGxfv1664hsXhITEzVt2jQ9/vjjGj58uHx8fNS4cWPNnz9fmzdv1sqVK2+rLwBzycrK0kcffaRHH320WP1B0Lx5c23evFljxoxR+fLlXV2OQwYPHqxz584pMjLS1aUUeQTZYiwtLU0zZ85U/fr15e3trSpVqmjkyJG6cOGCTZ+XX35ZderUkZeXl8qUKaPHH39csbGxLqzc1uXLlzVq1ChVqlRJnp6eCgsL04QJE5SUlCRJWr58uSwWi86cOaMVK1ZYvwb6/fffc12nI/sm29KlS1WrVi15e3urSZMm+vnnn22Wf/fdd7JYLIqOjrZpT01NlcVi0YwZM6xtf533eqv15lb3wIED5efn57QR7XXr1ik5OVl9+vSxae/YsaP8/f1twml++haUv+7DJUuWqHbt2vLw8LDuv8OHD2vgwIEKCQmRp6enatWqpddff91utPjw4cMaPHiwKlasKB8fHzVq1EgffvihTb/o6Gh17NhRAQEB8vX1VdOmTfXFF1/YrCd7HmdKSopGjRql4OBgBQQEqF+/foqLi7PpGxMToyFDhqhy5cry9fVVvXr19I9//MP6Wh48eLBmzZolSSpRooTdazn7ua5du6annnpKZcqUUWhoqCSpf//+uueee+z216JFi2SxWHT06FGHt/9WdeRnP69Zs0aNGzeWt7e3atSooU8++STXY5sbR/ZtTnMv73R/S469BiRp9erVCg8Pt9nOnOYC53UMs49V9s3f31/333+/3Xv9r+sYNmyYgoKCFBISoldeeUXSn39wDhs2TMHBwSpVqpTGjh3r0Mje3r17dfr0abVv396m/Y8//sj163aLxaLz58/fct234uh+LuqcdWzCw8MVHBysVatWFfYmmA5BtphKT09Xp06d9O6772rq1KmKjY3VL7/8oiZNmmjJkiWSpIyMDD300ENasGCBZs+erUuXLik6OlqxsbFq06aN4uPjnVbPqlWrVLJkSXl7e6t+/fqaPXu2bty4ccvHJScn64EHHlBUVJSWLl2qy5cva8GCBVq6dKk6duyo9PR09e/fX4ZhqFKlSurXr5/1a6AGDRrc9r7J9u2332rfvn3atGmTjh49qjJlyqh37953vG9uZ72XLl1Su3bttGnTJm3evNlpc4z37dsnSapdu7ZNu7u7u2rWrKnffvvttvoWtC+//FK///67fvrpJ/3yyy8KCAjQb7/9pmbNmikuLk7r169XXFyc3nzzTb322msaO3as9bF79uxR06ZNdebMGa1atUqXL1/W559/rl27dlm3ccOGDWrXrp2Cg4O1d+9enTx5Ut27d9egQYP09ttv29Uzfvx4de7cWSdPnlRkZKR+/PFHm+eUpJ49e+rw4cP64YcfFBcXp2+//Vbe3t7WPwA+//xzTZ06VdKfr9PcXsujR49Wly5ddOTIEZs/lBx1q+2/VR2O7ud169apZ8+eatq0qY4cOaJNmzZp586d+f4jzJF9m5M73d+OvgaioqLUu3dvtWzZUkePHtXmzZu1d+/ePLczp2M4fPhwaw2ZmZk6ePCgHnjgAT3yyCN2fyhL0rhx4/TYY4/p9OnTmjt3rqZPn64lS5ZoxIgRevjhh3Xy5EktXLhQ8+fP1/z582+5v7Kfo2nTpjbtderUsfuaPfv/rdDQUAUEBNxy3XlxdD9nDxo4crvTk2Qff/xxeXp6KigoSB07dtQPP/yQr8c749g0a9ZMmzdvvqPtuCsYKJYWLFhgSDJWrVqVa5+PP/7YkGSsXbvWpv3ixYuGr6+vMXPmTMMwDOPEiROGJOOjjz6y9smpLTcDBw40fvzxRyM+Pt44ffq08frrrxteXl5Gx44djaysrDwf+/bbbxuSjJ9++smmfeXKlYYkY8mSJda2SpUqGf369btlPY7sm927dxuSjPbt29u079+/35BkLFy40Nq2evVqQ5KxZcsWm74pKSmGJGP69Om3td7svsuWLTP2799vVK1a1QgPDzdiYmJuuY03q1SpktGiRYscl40aNcqQZJw7d85uWfv27Y3g4ODb6ltQsvdLq1at7Ja1a9fOCA0NNZKTk23aFyxYYLi5uRnHjh0zDMMw2rRpY1SoUMG4du1ars/TpEkTIywszLhx44ZNe69evQw/Pz8jPj7eMAzDmDhxoiHJWLx4sU2/KVOmGO7u7kZiYqJhGIYRHx9vSDLmzZuX5/ZNnTrVkGSkp6fbLct+rvnz59st69evnxEWFmbX/tFHHxmSjCNHjljbHNn+vOpwdD83btzYqFu3rt37vHHjxoYk4+DBg7k+/1+391b71jAMY+bMmYYkIyUlxTAM5+xvR18DjRo1MurXr2+3nc2aNbPbzryOYW4aNGhgDB061G4dn3zyiU2/9u3bG35+fsaiRYts2jt27Gg0bNjwls/z3HPPGZJs9mtOrly5YtSqVcsICAgw9u3b5/B2GIb9cTIMx/dz9v+1jtyWLVuW4/Nv3LjRkGS89957OS7ftGmT8cILLxj79u0zrl27ZuzZs8fo2bOn3edNbpx5bIYOHWq4u7vf8nPybseIbDEVFRUlf39/9ejRI9c+q1evVkBAgDp16mTTXrZsWd17773atGmTU2pZunSpOnTooMDAQIWGhur555/XP/7xD/3444+3nP+zYcMGBQUFqV27djbtvXv3lru7uzZs2JDvehzZN9m6detmc79OnTry8PDQ8ePH8/28t7veDRs2qHXr1qpfv762bNmiypUr39Fz58ZisTjcnp++f3XzFSzGjRuXZ3teevbsaXM/KSlJmzZtUrdu3eTr62uzrEOHDsrKytKWLVuUkJCgrVu3qlevXvLz88tx3VevXtWuXbvUo0cPlShRwmbZo48+quTkZG3bts2m/eZj2qBBA2VmZurUqVOSpMDAQN1zzz168803tWTJkjv6Ovbmbc8PR7Y/L47u56tXr2r37t3q3r273euiV69e+XrOW+3bnNzp/nb0NXD16lXt3btX3bp1s9vOvP6PyekYpqam6uWXX1b9+vXl4+NjM9Xh5qkhkvTQQw/Z3K9Tp46Sk5Pt2uvWrevQ/1nZ3wiVLFky1z5paWnq3bu3jh8/rq+++koNGza85Xrzkp/3Wvfu3e1GhnO79e/f/7bqeeCBB/Taa6+pYcOG8vPzU6NGjfT111/r3nvv1YQJExz6JlFyzrEJCAhQZmamdSoMckaQLaYuXryoihUr5tnn/PnzSkxMlKenpzw8POTu7i43NzdZLBZt27ZNV65cKbD6sr8Wz+nrsr+6cuVKjpPzPTw8VKZMGV2+fDnfz+3IvslWoUIFm/tubm7y8/NzaGqBYRhOWe+3336rxMREjRo1Ks8PmNtVunRpSX9+oNwsPj5ewcHBt9W3oN18JYaLFy8qKytLH374ofX1nP2azr6SxZUrV3T58mXrVJTcZM+/zOm1l93219eem5ubypUrZ9Mv++vWvx7T77//Xvfdd59Gjx6tChUqqE6dOvrHP/6ha9euObzdFovF7vWTl5tfh45sf14c3c/Z/3+EhITYrSOnttw4um9zcif729HXQPZ23lxjbm1S7sdw+PDhevPNNzVjxgzFxMQoMzNThmGoZcuWdtcUzWm/+Pv7S7L//8Xf31/Jycm3nCcbFBQkSbkGJ8MwNGzYMG3ZskUffPCBOnbsmOf6HJHf95oreHh4qFu3brp69aoOHDhwy/7OOjaJiYlyd3e3PhY5I8gWU2XLltXZs2fz7FOmTBmVL19eGRkZysjIUGZmprKysqx/0e7atauQqs1dcHBwjidgZWZm6sqVKypTpky+1+nIvsl2qxFG6c+RH8n+P/+8TphzZL3Z3nrrLT388MPq06dPgZxMlT2icujQIZv2zMxMHTlyRPfee+9t9c3J+fPnbUZN5s6dm2d7Xm4evSldurQsFosmTJhgfT3f/JqeMGGCypQpI4vFkufxyQ7kOb32stv++tpz9HjWrl1b33zzjeLj4/Xrr7+qT58+evXVVzV8+HCHHi/9+SHp7u5u1x4YGJhjALl5Ox3Z/rw4up+z/+jJax86Ij/vlZvdyf529DWQ3e/ixYt2/XJqk3I+hjdu3NCKFSs0YsQIPfbYYypTpozc3P78iD5x4oTdOvLaL7e7z8LCwiQp19HradOmadmyZXrxxRfz9ZrNS37ea4U5R/ZOOOvYnDt3TlWqVLmj98DdgCBbTHXv3l1JSUlavXp1rn169Oih8+fPu2QyefaUgvvvvz/PfhEREbp69ardNIfIyEhlZGQoIiIi38/tyL7Jj2rVqkmS3VUSvvvuO6es39PTUytWrNDgwYPVv39/p/+aWufOneXn56dvvvnGpv3HH39UUlKSzbUT89O3sAUFBalNmzZatWpVnr8aFBgYqDZt2igyMlLJyck59ilVqpQaN26s7777zm6k5Ouvv5avr+8tf2AiL56enmrZsqVmz56tLl262LwHs7/uz++vKlWvXl1xcXF2f6StWbPG5r4j259XHY7u51KlSik8PFxr1qyxGxUu7DOxb2d/O/oaCA4OVqNGjRQVFWW3nbfzf8DNP0Lwww8/5Cv434k2bdpIknbu3Gm37OOPP9arr76qAQMGWM/Ad4aCfq85Q0ZGhqKiohQUFKR69eoV2vPu2LFDbdu2LbTnMyuCbDE1bNgwtW3bVk8++aSWLVumK1euKDY2Vh999JHmzJlj7dO+fXsNGDBAy5Yt08WLF5WYmKidO3dqwoQJDp3leiuffPKJJk6cqF27dunatWuKjY3VO++8o3/+85/q0KHDLefKPfXUU6pdu7aGDh2qzZs3KykpSevWrdPo0aPVtGlTDRw4MN81ObJv8qNSpUrq2rWr5s6dq+3btysxMVHLly/XwYMH872u3Li7u+vjjz/W+PHjNWLECL355ptOW3dgYKBefvllffbZZ1q0aJFSUlK0e/dujRkzRvfff7/69u17W31d4f3339eFCxfUq1cv7dixQ8nJyYqNjdV3332nLl26WEch582bp6SkJHXr1k07d+5UcnKyfv/9d40aNUp79uyRJL322ms6ffq0nnjiCZ08eVKXLl3SP//5T3377beaMWOGdSTeUQcOHFDPnj31/fff68KFC0pJSdHGjRu1bds2mzng2WfMr169Ol8XQ3/88cfl7e2tsWPH6sKFCzpz5ozGjBlj/UPrrxzZ/rzqcHQ/z5o1S3/88YdGjhyp2NhYnT17VmPHjrVebqogOWN/O/oaePXVV3XgwAGNHj1asbGxOnfunJ577rl8TQHx9PRU586d9fHHH2vr1q26du2avv/+e82YMUP33XefE/bIrd17772qUqWKfvrpJ5v2HTt2aOTIkWrTpo2WLFni9BFCR/dzYcyRHTJkiD7++GMdPXpUKSkp+u2339S3b1/t3btXb731ljw9PZ256bnas2eP4uLi7mg+/F2j4M4jg6ulpKQY06dPN2rVqmV4enoaVapUMUaOHGlcuHDB2ufGjRvG66+/bjRq1Mjw9vY2goKCjBYtWhhz5841kpKSDMO4s6sWJCcnG/PnzzdatGhhBAQEGD4+PkZ4eLjx2muvGWlpaQ5tx8WLF42nnnrKqFChguHh4WGEhoYazz77rJGQkGDTz9GrFjiyb/56xYCbBQYGGn//+99t2s6ePWs9wzY4ONh45plnjISEhFyvWuDIenPrO2vWLEOS8eKLL+a5jUOGDMn1jN6pU6fa9f/oo4+M+vXrG56enkb58uWN0aNHW88WvpO+zpbXPjQMwzh+/Ljx5JNPGpUrVzZKlChhhIaGGr169TLWrVtn0+/gwYNGv379jLJlyxre3t5GeHi4sXDhQiMzM9PaZ9OmTUb79u2NkiVLGt7e3kaTJk2Mf//73zbrmThxouHu7m5Xx81Xs8jKyjJWr15tdO3a1QgJCTH8/PyMunXrGjNmzDCuX79ufVxWVpYxZswYo1y5cobFYjEkGb/99luez5VtzZo1RoMGDQxPT0+jdu3axsqVK3O8aoEj259XHfnZz6tWrTIaNWpkeHp6GtWqVTMWLVpkLFmyxOGrFjiybw3D/mx4Z+xvw3DsNWAYhhEZGWndzurVqxtLliyxXiHl+PHjt9wmw/jz/7pBgwYZZcqUMfz9/Y1u3boZx48fNx588EGbq47kto5JkyYZOX2s53Vlhpu98sorRmBgoM1VBT777LM8rxCQ01VMcpPTVQsMw/H9fLvc3d1zrf+vr8OjR48azzzzjPWzITg42HjooYeMH374waHncdaxmThxolGhQgW7KznAnsUw8jgjBQAA3JaXXnpJr776qhITE2/rChGuEBcXpxo1amjOnDl66qmnXF3OXenatWuqWrWqXnzxRY0fP97V5RR5BFkAAJzMMAw1bNhQvr6++u9//+vqcvJl7ty5evvtt3XkyBG7ObsoeK+++qo++eQT7d+/3+6kVtgjyAIAcAfS0tI0dOhQjR8/XvXr11dsbKxmzpypL774QqtXr1bXrl1dXSJQbHGyFwAAd8DLy0t9+vTR+PHjVbFiRTVp0kQnT57Ud999d1eE2PXr19/yclhjxoxxdZkophiRBQAAgCkxIgsAAABTIsgCAADAlDxcXUBRkJWVpbNnz8rf35+fggMAACgghmEoKSlJFStWtP4M850gyEo6e/ZsofzSDAAAAKSYmBhVrlz5jtdDkJXk7+8v6c+dGhAQ4OJqAAAAiqfExESFhoZas9edIshK1ukEAQEBBFkAAIAC5qypnJzsBQAAAFMiyAIAAMCUCLIAAAAwJYIsAAAATIkgCwAAAFMiyAIAAMCUCLIAAAAwJYIsAAAATIkgCwAAAFPil72Am3Seucahfute6lbAlQAAgLwwIgsAAABTIsgCAADAlAiyAAAAMCWCLAAAAEyJIAsAAABTIsgCAADAlAiyAAAAMCWCLAAAAEyJIAsAAABTIsgCAADAlAiyAAAAMCWCLAAAAEyJIAsAAABTIsgCAADAlAiyAAAAMCWCLAAAAEzJw9UFZGVl6cKFCypXrpzc3d1tlhmGoWPHjtk9JiQkRP7+/nbt586dk5eXl4KDgwusXgAAABQNLhuRPX36tIYPH66goCA1adJEJUuW1KhRo5SWlmbtk5ycrJo1a6pdu3bq0qWL9fbjjz/arGv37t2qV6+eatWqpQoVKqhTp066dOlSYW8SAAAACpHLguzWrVvVqlUrnT9/XufOndPevXsVGRmpKVOm2PVduXKljh49ar09/PDD1mUpKSnq2bOnWrduratXr+rixYu6evWqhg4dWohbAwAAgMLmsqkFAwYMsLlfq1YtDRgwQOvWrbPre/36dcXGxqpixYqyWCw2y1avXq2zZ8/q1VdflYeHhwIDAzVt2jT17t1bp0+fVpUqVQp0OwAAAOAaRepkrz/++EOVKlWya3/ooYcUHh4uPz8/Pfvss7p+/bp12Y4dO1S9enWVK1fO2tamTRtJ0s6dOwu+aAAAALiEy0/2yrZy5UqtXbtW33//vbXNzc1Nr7/+usaMGSMfHx/t3LlTPXr0UHp6uv71r39Jki5fvqzSpUvbrKtUqVJyc3PT5cuXc3yutLQ0m7m4iYmJBbBFAAAAKEhFYkR248aNeuKJJzR79mx16tTJ2u7r66vnn39ePj4+kqSmTZtq8uTJ+uSTT5SVlSXpz7Cbnp5us77MzExlZWXZXQUh2+zZsxUYGGi9hYaGFtCWAQAAoKC4PMhu2rRJPXr00NSpUzVp0qRb9r/nnnuUkpKiixcvSpJCQ0N1/vx5mz7Z9ytXrpzjOqZMmaKEhATrLSYm5g63AgAAAIXNpUF2y5Yt6tatmyZNmqRp06bZLc/IyLBri46OVmBgoMqWLStJ+tvf/qbY2Fjt37/f2icqKkpeXl5q2bJljs/r5eWlgIAAmxsAAADMxWVzZLdt26auXbuqb9++GjBggI4ePSpJcnd3V9WqVSVJ7777rk6fPq2ePXsqKChIUVFRmjdvnubMmWOdNvC3v/1N7dq10+OPP6533nlHcXFxevHFFzV+/HgFBga6avMAAABQwFwWZKOjoxUSEqLNmzdr8+bN1vaAgADt2rVLkvTss89q8eLFmjNnji5duqRq1appzZo16tixo826vvnmG82YMUNjx46Vl5eXNcgCAACg+LIYhmG4ughXS0xMVGBgoBISEphmAHWeucahfute6lbAlQAAULw4O3O5/GQvAAAA4HYQZAEAAGBKBFkAAACYEkEWAAAApkSQBQAAgCkRZAEAAGBKBFkAAACYEkEWAAAApuSyX/YC7hb8wAIAAAWDEVkAAACYEiOywG1ydKQVAAAUDEZkAQAAYEoEWQAAAJgSQRYAAACmRJAFAACAKRFkAQAAYEoEWQAAAJgSQRYAAACmxHVkcdfguq8AABQvjMgCAADAlAiyAAAAMCWCLAAAAEyJIAsAAABTIsgCAADAlAiyAAAAMCWCLAAAAEyJIAsAAABTIsgCAADAlAiyAAAAMCWCLAAAAEyJIAsAAABTIsgCAADAlAiyAAAAMCWXBdnr169r7ty5atasmUqXLq3w8HDNnz9fhmHY9Dt79qz69u2r0qVLq1KlSpowYYLS0tLy3QcAAADFi4ernvjjjz/WqVOn9MEHH6hq1arasmWLBg8erOvXr+v555+XJGVmZqpbt24qXbq0tm/frqtXr+qRRx5RcnKyFi5c6HAfAAAAFD8W4+YhUBcaM2aMtm7dqt27d0uS1q5dq65du+rYsWOqVq2aJOnzzz/X0KFDdf78eZUpU8ahPreSmJiowMBAJSQkKCAgoOA2EC7VeeYaV5eQp3UvdXN1CQAAFChnZ64iNUc2Pj5eJUuWtN6Pjo5WWFiYNaBKUkREhDIzM7Vt2zaH+wAAAKD4cdnUgptt27ZNK1as0IcffmhtO3funMqVK2fTr2zZsrJYLDp//rzDfW6WlpZmM4c2MTHRWZsBAACAQlIkRmQPHz6s3r17a/DgwRo2bJjNMjc32xItFosk2ZwU5kifv5o9e7YCAwOtt9DQ0DveBgAAABQul4/IHj16VO3bt1eHDh308ccf2ywLCQnRpk2bbNouX74swzAUEhLicJ+bTZkyRRMmTLDeT0xMJMzC5fIzh5f5tAAAuDjIHjt2TO3atdMDDzygTz/91G5ktVWrVpozZ45iYmKsQXPjxo1yc3NT8+bNHe5zMy8vL3l5eRXglqEwFfWTuAAAQMFw2dSCEydOWEPsZ599Jnd3d7s+Dz30kOrUqaMxY8bo8uXLOnLkiKZPn64BAwaofPnyDvcBAABA8eOyIPvhhx8qJiZGX331lfz8/OTt7S1vb2+b8FmiRAmtWbNG169fV8WKFRUeHq7WrVtrwYIF+eoDAACA4sdl15HNyMhQRkaGXbvFYsnxa/+srCy7qQe30ycnXEfW3O7GqQXMkQUAmJGzM5fL5sh6eHjIw8Pxp3ckoN5OiAUAAIA5kfwAAABgSgRZAAAAmBJBFgAAAKZEkAUAAIApEWQBAABgSgRZAAAAmBJBFgAAAKZEkAUAAIApEWQBAABgSgRZAAAAmBJBFgAAAKZEkAUAAIApEWQBAABgSgRZAAAAmBJBFgAAAKZEkAUAAIApEWQBAABgSgRZAAAAmBJBFgAAAKZEkAUAAIApEWQBAABgSvkOsrGxsTp58mS+lwEAAADOlO8gu2zZMr3//vv5XgYAAAA4k1OnFsTHxysgIMCZqwQAAABy5OFox6ioKH3xxRc6ePCgUlJSdP78eZvlycnJ+vHHH/Wf//zH6UUCAAAAN3M4yLq5ucnDw0Nubm7Wf/9VaGioPv30U3Xq1MnpRQIAAAA3czjIdunSRV26dNGmTZuUlJSk7t27F2RdAAAAQJ4cDrLZHnzwwYKoAwAAAMiXfAdZSdq9e7fef/99nThxQjdu3LBZNmDAAI0ePdopxQEAAAC5yXeQPXXqlNq2bau2bduqVatWKlGihM3ymjVrOq04AAAAIDf5DrLr169X27ZttXbt2oKoBwAAAHBIvq8j6+vrq9DQ0IKoBQAAAHBYvoPsAw88oM2bN+vKlSt3/ORJSUlasGCB2rVrp2HDhtktT0lJUcuWLe1ua9assemXmpqqV199VREREerWrZv+/e9/33FtAAAAKNryPbVg//79slgsqlOnjtq3by9/f3+b5Z07d9Zjjz12y/XcuHFDtWvXVvfu3eXv76/ffvvNrk9mZqa2b9+uRYsWqX79+tb26tWr2/Tr37+/Dh48qFmzZunq1asaPXq0zp49q8mTJ+d38wAAAGAS+Q6yN27cUKNGjaz3r127ZrfcESVKlNChQ4fk7++vcePG6ezZs7n2rV+/vlq2bJnjsl9//VWRkZHasWOHmjZtKunPXxmbNm2annnmGfn5+TlUDwAAAMwl30G2e/fuTvkxBIvFYjeam5tJkybJzc1N1atX14gRI9SiRQvrsg0bNqh8+fLWECtJvXr10vjx47Vt2zZFRETcca0AAAAoevI9R7awNW3aVMOGDdPkyZNVsmRJtWnTRl9++aV1+alTp1SxYkWbx1SqVMm6LCdpaWlKTEy0uQEAAMBc8j0iGxkZqSVLluS6vHfv3ho6dOid1GTl6+urX375xXqt2s6dOysjI0MTJkxQ3759JUnp6eny8vKyeVyJEiXk5uam9PT0HNc7e/Zsvfzyy06pEQAAAK6R7xHZkiVLqnLlyja3UqVKadeuXdq7d6+CgoKcV5ybm90PLnTo0EGxsbHWqyYEBwcrLi7Opk98fLyysrJUunTpHNc7ZcoUJSQkWG8xMTFOqxkAAACFI98jshERETnOO71+/bpatWqlatWqOaWw3Jw/f15ubm7y9vaWJDVp0kTvvfeerl69qlKlSkmStm/fLklq3Lhxjuvw8vKyG8UFAACAuThtjqyvr6969eql9evXO2uVWrVqlXbs2GG9f+TIEc2ZM0c9evSwXo2gV69eCgoK0qxZsyT9edWEOXPm6MEHH7S7TBcAAACKj3yPyObl+PHjatCggcP9Bw0apGPHjunUqVNKSkqyXmJr48aN8vHxUfXq1fX000/r6NGjCgwM1LFjx9SvXz/NmzfPug5/f399/fXX6t+/v1asWKHk5GSFhYUpMjLSmZsGAACAIsZiGIaRnwds2bLF7pe1MjMz9dtvv+nnn3/W//73P5sfL8jLvn37dP36dbv25s2by83t/w8WX7x4UZcuXVLVqlXl6+ub47oyMjJ08OBBeXl5qVatWvnYIikxMVGBgYFKSEhQQEBAvh6LgtN55ppbd0Ke1r3UzdUlAABg5ezMle8R2fPnz2vnzp22K/HwUJUqVRQdHe1wiJWke++916F+5cqVU7ly5fLs4+HhoYYNGzr83AAAADC3fAfZxx57zKGfoAUAAAAKUpH/QQQAAAAgJ7cVZNPT0/Xee++pU6dOqlevniIiIjR79mylpKQ4uz4AAAAgR/kOsoZhqHPnzpo+fbqqVaumQYMGqV69enr33XfVqlUrpaWlFUSdAAAAgI18z5H98ccfdeTIER08eFAhISHW9lmzZqlVq1ZasWKFnnjiCacWCQAAANws3yOyf/zxh7p06WITYiUpICBADz/8sA4dOuS04gAAAIDc5DvIlitXTvv27VNWVpbdsl27dqls2bJOKQwAAADIS76DbLdu3RQTE6MePXpo1apV+t///qeoqCj17dtX0dHR6tu3b0HUCQAAANjId5D19/fXzz//LMMw9PDDD6tp06bq0aOHLly4oJ9//lkVK1YsiDoBAAAAG/k+2UuSatWqpaioKKWnp+vs2bMqX768vLy8nF0bAAAAkKvbCrLZSpQoobCwMGfVAgAAADjstq4j26dPH23fvt2m/ciRI+rYsaPS09OdVhwAAACQm3wH2fXr1+v69etq0aKFTXvNmjVVoUIFrVy50mnFAQAAALnJd5A9dOhQrtMJwsLCdPDgwTsuCgAAALiVfAfZatWqaePGjUpNTbVpz8zM1Lp165gzCwAAgEKR75O9OnXqpBIlSqhDhw4aP368wsLCdPbsWc2fP1+xsbFcRxYAAACFIt8jsh4eHlq3bp0CAgLUt29fNWvWTL1791Zqaqo2bNiggICAgqgTAAAAsHFbl98KDQ1VVFSUkpKSdPHiRZUuXVpBQUFOLg0AAADI3R1dR9bf31/+/v7OqgUAAABwWL6nFgAAAABFAUEWAAAApkSQBQAAgCkRZAEAAGBKBFkAAACYEkEWAAAApkSQBQAAgCkRZAEAAGBKBFkAAACYEkEWAAAApkSQBQAAgCkRZAEAAGBKBFkAAACYkoerC9i2bZtWrFihcuXKacqUKTn2iYqK0oYNG+Tt7a1HH31UjRs3vq0+AAAAKD5cNiKbkZGhxo0ba+zYsdq2bZu+/vrrHPs9//zzGjx4sHx8fBQfH68WLVpo5cqV+e4DAACA4sVlI7IWi0VLlixReHi4xo0bp+joaLs+hw4d0ltvvaVVq1ape/fukiR/f38988wz6tOnjzw8PBzqAwAAgOLHZSOy7u7uCg8Pz7PPd999p6CgID300EPWtsGDB+vChQv673//63AfAAAAFD9F+mSvI0eOKDQ0VO7u7ta2qlWrSpKOHj3qcJ+bpaWlKTEx0eYGAAAAcynSQTYlJUX+/v42bX5+fnJ3d9f169cd7nOz2bNnKzAw0HoLDQ0tmA0AAABAgSnSQTYgIEDx8fE2bYmJicrMzFRAQIDDfW42ZcoUJSQkWG8xMTEFUT4AAAAKUJEOsvXr19eJEyeUmppqbTtw4IB1maN9bubl5aWAgACbGwAAAMylSAfZXr16KSsrS0uWLLG2vf/++6pbt64aNWrkcB8AAAAUPy69NtXMmTMVGxurX375RefPn9fIkSMlSfPmzZOXl5cqVKig+fPna8yYMYqKilJcXJz++OMPRUVFWdfhSB8AAAAUPy4NsrVr11bZsmXtLsP11ysQPPnkk2rfvr22bNkiLy8vdezYUaVKlbLp70gfAAAAFC8WwzAMVxfhaomJiQoMDFRCQgLzZYuQzjPXuLqEu8a6l7q5ugQAwF3A2ZmrSM+RBQAAAHJDkAUAAIApEWQBAABgSgRZAAAAmBJBFgAAAKZEkAUAAIApEWQBAABgSgRZAAAAmBJBFgAAAKZEkAUAAIApEWQBAABgSgRZAAAAmBJBFgAAAKZEkAUAAIApEWQBAABgSgRZAAAAmBJBFgAAAKZEkAUAAIApEWQBAABgSgRZAAAAmBJBFgAAAKZEkAUAAIApEWQBAABgSgRZAAAAmJKHqwvA3afzzDWuLgEAABQDjMgCAADAlAiyAAAAMCWCLAAAAEyJIAsAAABTIsgCAADAlAiyAAAAMCWCLAAAAEypSF9HNiMjQ4sWLbJrb9eunWrXrm3TdurUKW3evFne3t7q0KGDSpUqVVhlAgAAwAWK9IhsamqqRo0apR9++EF79uyx3uLi4mz6LV68WHXr1tWKFSs0d+5c1ahRQ9u3b3dR1QAAACgMRXpENtsLL7ygli1b5rjs3LlzGj16tN5++22NGjVKkjR48GANGzZMBw4cKMwyAQAAUIiK9Ihsts2bN+vTTz9VdHS0MjMzbZZFRkbKzc1NQ4cOtbaNGTNGBw8e1N69ewu5UgAAABSWIj8i6+7urqioKFWoUEFbtmxRmTJl9M0336hq1aqSpP3796tq1ary8fGxPqZevXrWZY0aNbJbZ1pamtLS0qz3ExMTC3grgKKt88w1DvVb91K3Aq4EAADHFekRWU9PT23btk0///yzli1bpkOHDsnT01PDhw+39klMTFRQUJDN4wICAuTu7p5rQJ09e7YCAwOtt9DQ0ILcDAAAABSAIj0i6+npqaZNm1rv+/n56ZlnntGwYcOUmpoqb29v+fj4KCkpyeZxycnJyszMlK+vb47rnTJliiZMmGC9n5iYSJh1AkdH9QAAAJyhSI/I5sTb21uZmZnW0daaNWsqJibGZu7siRMnJEk1atTIcR1eXl4KCAiwuQEAAMBcinSQPXnypLKysmzavvjiC1WvXl3lypWTJHXv3l3x8fFau3attc/nn3+ukJAQNW/evFDrBQAAQOEp0lMLfvnlFz3yyCPq3LmzgoKCFBUVpX379mnlypXWPrVr19bEiRP1xBNPaPTo0YqLi9OHH36opUuXysOjSG8eAAAA7kCRHpEdOHCgli9frqCgIF26dEn9+vXT0aNHFRERYdPvjTfe0Geffabr168rMDBQ27dvV9++fV1UNQAAAAqDxTAMw9VFuFpiYqICAwOVkJDAfNk7wMlexR+X3wIA3AlnZ64iPSILAAAA5IYgCwAAAFMiyAIAAMCUCLIAAAAwJYIsAAAATIkgCwAAAFMiyAIAAMCUCLIAAAAwJYIsAAAATMnD1QUAMA9Hf72NXwADABQGRmQBAABgSgRZAAAAmBJBFgAAAKZEkAUAAIApEWQBAABgSly1ALfk6JnqAAAAhYkRWQAAAJgSQRYAAACmRJAFAACAKRFkAQAAYEoEWQAAAJgSQRYAAACmRJAFAACAKXEdWQBO5+i1h9e91K2AKwEAFGeMyAIAAMCUCLIAAAAwJYIsAAAATIkgCwAAAFMiyAIAAMCUCLIAAAAwJS6/dRdz9BJJAAAARREjsgAAADClYjMim56ert9//13e3t6qW7euq8txKUZaAQDA3aBYBNmNGzdqwIAB8vHxUVJSkipVqqRVq1YpLCzM1aUBAACggFgMwzBcXcSdSExMVLVq1fR///d/eu2115Senq7OnTsrIyNDmzdvdngdgYGBSkhIUEBAQAFXXPAYkUVxw0/ZAkDx4OzMZfo5spGRkUpMTNSUKVMkSSVKlNCkSZO0ZcsWHT161MXVAQAAoKCYfmrB7t27Va1aNQUFBVnbmjdvbl1Wo0YNu8ekpaUpLS3Nej8hIUHSn38lFAcZqdddXQLgVBFTVzrU75tJnQu4EgDAncjOWs6aEGD6IBsXF6fSpUvbtAUFBcnNzU1xcXE5Pmb27Nl6+eWX7dpDQ0MLpEYAhSPwVVdXAABwxJUrVxQYGHjH6zF9kC1RooRSU1Nt2m7cuKGsrCx5enrm+JgpU6ZowoQJ1vvx8fEKCwvT6dOnnbJTcecSExMVGhqqmJiYYjFvuTjgmBQ9HJOiheNR9HBMip6EhARVqVJFwcHBTlmf6YNsWFiYVq9ebdMWGxsrSapSpUqOj/Hy8pKXl5dde2BgIC/0IiYgIIBjUsRwTIoejknRwvEoejgmRY+bm3NO0zL9yV4dO3bUhQsX9N///tfaFhkZqZIlS6pVq1YurAwAAAAFyfQjsi1atFCfPn00aNAgvfLKK4qLi9NLL72k6dOny9fX19XlAQAAoICYPshK0rJlyzR37lwtXrxYXl5eWrhwoQYNGuTw4728vDR9+vQcpxvANTgmRQ/HpOjhmBQtHI+ih2NS9Dj7mJj+BxEAAABwdzL9HFkAAADcnQiyAAAAMCWCLAAAAEypWJzslR8pKSnavXu3KlWqpLCwMJtlcXFxOnDggN1jmjdvnuuPK+DOZGVl6ejRo8rKylK1atVy3c/nzp3TmTNnVL16daddRBk5S09P16FDh+Tn56ewsDC7a/0dOHDA7lfzSpUqpfr16xdmmXeVGzdu6NChQ/Lx8VHVqlXl7u5u1ycrK0sHDhxQZmamGjRokGMfOE9ycrIOHz6s0qVL212zPD09Xdu3b7d7TN26de1+iRLO9+uvv8rT01P33Xef3bKEhAQdOXJE5cuXV+XKlV1Q3d0nNTVVO3fuVJkyZVSnTh1re0xMjE6dOmXT18PDQy1btszfExh3iQsXLhjjxo0zKlSoYHh7exuTJk2y6/PNN98Ybm5uxv33329zu3jxogsqLv7mzZtnVK5c2ahZs6ZRo0YNo2zZssbSpUtt+mRkZBhDhw41vL29jXr16hleXl7Gyy+/7KKKi7eUlBTjhRdeMIKDg40GDRoYFStWNGrVqmVs3brVpl+3bt2MihUr2rxHJk6c6KKqi7fMzExj2rRpRtmyZY3w8HCjYsWKRpUqVYzvv//ept/+/fuNGjVqGOXLlzcqVapkVKlSxdixY4eLqi7e4uLijOHDhxvBwcFGkyZNjODgYCM8PNw4cOCAtc+5c+cMSUbjxo1t3icbN250XeF3ifnz5xtubm5G9erV7Za98847hre3t1G3bl3Dx8fHePTRR43U1FQXVHl3GTFihOHm5mY88sgjNu0zZ840/P39bd4jXbp0yff675oR2VOnTik0NFS///67HnjggVz7+fj4KDo6uhAru3vFxcVpx44dKl++vCTpvffe05AhQ3Tfffepdu3akqS5c+dq1apV+v3331W9enVt2rRJERERatq0qbp27erK8oudhIQElS5dWqdOnVLJkiWVmZmpp59+Wr169dKZM2dsLpXy2GOPae7cua4r9i6RlpamUqVK6dSpU/Lx8ZFhGBo/frz69u2ruLg4ubu7KysrS3379lV4eLi+/PJLWSwWDR06VI8++qgOHz7Mt0lOdv78eUVERGjhwoVyc3NTWlqaevTooREjRth9dvz73/9WgwYNXFTp3Wffvn2aPXu2hg4dqk2bNtks27p1qyZMmKDvvvtOXbt21ZkzZ9SsWTPNmjVL//znP11UcfG3cuVK7dy5UxERETkur1ev3h1nrrtmjmyzZs00YcIEh76WPnjwoH777TelpqYWQmV3rxkzZlhDrCSNHDlSmZmZNl/JLV68WAMHDlT16tUlSQ8++KDatm2rxYsXF3q9xV1ISIheeOEFlSxZUpLk7u6up556SpcvX9ahQ4ds+l67dk07d+7U6dOnZXAFvwLj4+OjCRMmyMfHR5JksVj04IMPKjExUUlJSZKk7du3a//+/Zo2bZosFosk6aWXXtKpU6f0008/uaz24qpu3brq37+/dcqNl5eXWrVqZf1p9L86ffq0du3apYSEhMIu865z/fp19e/fX++9954qVKhgt3zx4sW67777rAMglStX1tChQ/ksKUAnTpzQs88+q6VLl+b6B3VGRob27t2rQ4cOKSMj47ae564Jso5KTk5Wjx491KdPH5UqVUozZ850dUl3jf/9738yDEM1atSQ9Oe8moMHD9rNc2revLl2797tihLvOjt27JC7u7vuuecem/bPP/9cw4cPV3h4uOrWratff/3VNQXeJY4dO6YtW7boiy++0IsvvqgXXnhBQUFBkqTdu3fLw8ND9957r7V/9lxy3icFZ9++fdq0aZMWLFighQsX6uWXX7br8+STT+qJJ55Q2bJlNWTIEF27ds0Fld4dnnnmGbVp00a9e/fOcfnu3btz/CyJjY3VpUuXCqHCu0tGRoYGDBigadOmqW7durn22717twYOHKiIiAiVL19en332Wb6f666ZWuCI0NBQ7dy50/piX7t2rXr27KnKlStr2LBhLq6ueLt27ZpGjBihiIgItW7dWpIUHx8vwzDsTo4oXbq03clGcL5jx45p6tSpevbZZxUQEGBtHzx4sL744gsFBAQoLS1NTz31lHr37q0DBw5wIksB+eabb/Sf//xHJ06cUIUKFdSvXz/rsri4OAUHB1tHY7PxPilY8+fP1549e3T48GG1bt3a5qtTLy8vffnll3rsscckSYcOHVL79u01ceJELVy40FUlF1vLly9XdHS0du3alWufuLi4HD9LspeVLVu2QGu820ydOlVlypTR008/nWufZs2a6fjx49YT7+fOnauhQ4eqZs2a+TrhixHZv7jvvvts/mJ76KGH1KtXLy1fvtyFVRV/KSkp6tmzpwzDsNnXJUqUkCS7KR4pKSnM+ytgZ8+eVefOndW6dWvNmTPHZln//v2twdbLy0vz5s3TxYsX+Rq7AD333HP65ZdfFBsbqy5duqhdu3bWUaQSJUrkOA2K90nBWrhwobZv367Y2Fh5eXmpU6dOysrKkvTnVTyyQ6wk1a5dW+PHj9eKFStcVW6xFR8fr7///e8aNWqUdu/erejoaMXExCg1NVXR0dHWP+Zyep+kpKRIEu8TJ9u5c6fmzZunJ598UtHR0YqOjtbVq1d15coVRUdHKy0tTZLUuXNnm6tHjRs3TjVr1tTKlSvz9XwE2VsICQnJce4TnCM1NVU9e/bUhQsX9NNPP6lMmTLWZcHBwfL397fb/7GxsXaXu4HznD17Vu3atVPt2rX19ddfW/+gyE1gYKC8vb15nxQCNzc3TZw4UYmJidq6daskKSwszGbOrPTn5bouXbrE+6QQ+Pr6avTo0Tpw4IBOnz6da7+QkBAlJCQwvcDJUlNT1bBhQ3311VeaPHmyJk+erI0bN+rKlSuaPHmy9u/fL+nP90lOnyXu7u6qWLGiK0ovtlJSUtS0aVO9/fbb1mOSfe7R5MmTdfXq1VwfezuZi6kFf5GcnCw/Pz/r/YyMDP38889q2LChC6sqvrJDbGxsrDZu3Khy5crZLLdYLIqIiNCqVas0btw4SX9enzEqKkpPPPGECyou/s6dO6d27dqpevXq+s9//mM3UnHjxg25ubnJw+P//9exefNmpaamcnZ2Abj5/yRJOnr0qKT//7Vou3bt5OHhodWrV2vgwIGSpHXr1unGjRvq0KFD4RZ8F8jtmLi7u1vnLefU54cfflBYWJj1ZEo4R/ny5e3Oep82bZp1ukG2jh076vXXX1dqaqq8vb0lSZGRkXrggQdsrsiCO9e2bVu7Y9K9e3d5e3vrq6++srbd/D65dOmS9uzZo44dO+br+e6aIPvXC1Rfv35dZ86cUXR0tAICAqwnSTz55JMKCwtTmzZtlJ6ergULFujs2bP68ssvXVl6sfXwww9r+/btWrx4sY4cOaIjR45Iku655x7rhaqnT5+u1q1ba8yYMercubMWL16szMxMa7CF8yQkJCgiIkKZmZmaOHGiduzYYV3WoEEDBQUF6eLFi+rRo4eGDx+umjVr6tChQ3rllVfUtWvXXC+vgtu3fv16LViwQP369VPlypV1+PBhvfbaa+rYsaPatGkj6c8RjPHjx2vs2LFKS0tTiRIl9MILL2jEiBHWq33Aed555x0dO3ZMnTt3VnBwsHbu3KnXXntN48aNswbZd999VwcOHFC3bt3k7++vyMhILV++XEuXLnVt8XexkSNHasGCBerTp49GjRqlLVu2aO3atdq4caOrS7trtWvXTr1791bjxo11+fJlvf766ypXrlye82pzYjHukmvnxMXFqWfPnnbt9evXt06+T0tL00cffaSffvpJWVlZatCggZ599lkmgReQ7A/im40aNUqDBg2y3t+9e7fmzp2rM2fOqHbt2po0aZLdr7Lhzh09elRDhw7Ncdlbb72lFi1aSPrzJLD58+dr//79CgkJUadOnTRo0CC7k43gHFu3btWnn35qPdGrc+fO6t+/v80vd2VlZWnRokWKjIxUVlaWunbtqlGjRtmMnMN5vv76a33zzTfW6Rt9+/a1G0XK7nPlyhXVrFlTI0eOVL169VxU8d3lo48+0oYNG+zObzl37pzmzJmj33//XeXLl9eYMWPUqlUrF1V5d5k0aZI8PT1trgR19epVvf/++9q+fbt8fX3VvHlzPf300/L19c3Xuu+aIAsAAIDihZO9AAAAYEoEWQAAAJgSQRYAAACmRJAFAACAKRFkAQAAYEoEWQAAAJgSQRYAAACmRJAFABNZv369YmJi7mgdmzdv1vHjx51UEQC4DkEWAExi3759evzxx60/hXq7YmJiNHDgQPF7OADMjiALACbx4osv6umnn5a/v/8drWfAgAE6f/68vv32W+cUBgAuwk/UAoALJCUlac2aNXbtQUFB6tKli137iRMnVKNGDZ04cUJVqlSRJP3666/y8PBQnTp1tGfPHiUlJenBBx+Un5+fEhIS9Msvv8jb21utW7eWl5eXzfqmTp2q7du3a/369QWzgQBQCDxcXQAA3I2SkpLsRkR//vlnhYSE5Bhk16xZo2rVqllDrCTNmzdPf/zxh65evar69evr8OHDSk1N1axZszRjxgzVq1dPBw4cUHBwsH799Vd5enpaH9u+fXu9/vrrSkpKuuMRXgBwFYIsALhAxYoVtXz5cuv9qKgoff3111q4cGGO/Xft2qV69erZtR87dkz79u1T1apVlZqaqqpVq2rcuHHat2+fQkNDlZSUpLCwMH311VcaOHCg9XENGzZURkaG9uzZo7Zt2zp/AwGgEBBkAcDFDh06pIEDB2r69Onq1atXjn0uX76s4OBgu/ZOnTqpatWqkiRvb2+Fh4crODhYoaGhkiR/f3/raO1flSpVyrpeADArTvYCABdKSEhQz5491bFjR02dOjXXfiVLllRycrJde3Ygzebl5ZVjW2pqqk1b9rqYVgDAzAiyAOAiWVlZGjBggHx8fPTpp5/KYrHk2rdWrVo6ceKE0547e121atVy2joBoLARZAHARSZPnqydO3cqMjJSvr6+efaNiIjQ3r17de3aNac899atW1WrVi2bk8cAwGwIsgDgAlu3btUbb7yhvn376tdff9Xy5cu1fPlyff/99zn2b9OmjWrWrKmvvvrKKc+/YsUKDR8+3CnrAgBX4WQvAHABLy8v9evXT5cvX7a5DFdoaGiOl9+yWCyaNm2a3nzzTQ0ZMkQWi0WtW7eWt7e3Tb+2bdvanRTWrl07VatWzXp/9+7dOnLkiJ566innbhQAFDJ+EAEATGTcuHH6+9//rrp16972Oj744AOFhITo4YcfdmJlAFD4CLIAAAAwJebIAgAAwJQIsgAAADAlgiwAAABMiSALAAAAUyLIAgAAwJQIsgAAADAlgiwAAABMiSALAAAAUyLIAgAAwJQIsgAAADAlgiwAAABM6f8BlWX1u+9wj1QAAAAASUVORK5CYII=", - "text/plain": [ - "
" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], + "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "\n", + "from zagg.readers import chunk_z_range\n", + "\n", "n_bins, resolution = 128, 0.5\n", - "vec = tensor[0, 5] # the digest histogram for cell 5 of chunk 100\n", + "tensor = tensors[word_a]\n", + "vec = tensor[0, 5] # the digest histogram for cell 5 of chunk A\n", "\n", "# Reconstruct the window the reader used: floor of the chunk's trimmed minimum.\n", "# We recompute the floor here only to label the x-axis; the reader sets it internally.\n", - "from zagg.readers import chunk_z_range\n", - "\n", - "digests = [d for _, d in iter_csr_cells(read_csr(store, f\"{FIELD}/100\"))]\n", + "digests = [read_cell(store, FIELD, base_a + c) for c in (0, 5, 4095)] # chunk A's cells\n", "z_lo, n_bins_c, res_c = chunk_z_range(\n", " digests, n_bins=n_bins, resolution=resolution, bottom=0.05, top=0.95, fit=\"raise\"\n", ")\n", @@ -302,7 +312,7 @@ "ax.bar(centers, vec, width=res_c, align=\"center\", color=\"steelblue\", edgecolor=\"none\")\n", "ax.set_xlabel(\"z (m)\")\n", "ax.set_ylabel(\"count\")\n", - "ax.set_title(f\"cell 5 of chunk 100 -- reconstructed histogram (z_lo={z_lo:.0f} m)\")\n", + "ax.set_title(f\"cell 5 of chunk A -- reconstructed histogram (z_lo={z_lo:.0f} m)\")\n", "ax.set_xlim(z_lo, z_lo + 30) # zoom to the populated region\n", "plt.tight_layout()\n", "plt.show()" @@ -322,7 +332,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": null, "id": "22da8a78", "metadata": { "execution": { @@ -332,17 +342,7 @@ "shell.execute_reply": "2026-06-25T19:05:18.296964Z" } }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "dtype=uint16 -> tensor.dtype=uint16, peak bin count=295\n", - "dtype=uint32 -> tensor.dtype=uint32, peak bin count=295\n", - "dtype=float32 -> tensor.dtype=float32, peak bin count=295.212646484375\n" - ] - } - ], + "outputs": [], "source": [ "for dtype in (\"uint16\", \"uint32\", \"float32\"):\n", " t, _ = next(read_tensors(store, FIELD, dtype=dtype))\n", @@ -370,7 +370,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": null, "id": "fb8f0800", "metadata": { "execution": { @@ -380,18 +380,11 @@ "shell.execute_reply": "2026-06-25T19:05:18.323897Z" } }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "fit='raise' -> trimmed z-range [0, 400] (span 400) exceeds the fixed window 128 bins × 0.5 = 64.0; pass fit=\"degrade_resolution\" or fit=\"collapse_bins\" to adapt\n" - ] - } - ], + "outputs": [], "source": [ "wide = MemoryStore()\n", - "write_chunk(wide, FIELD, 1, {0: rng.uniform(0.0, 400.0, 5_000)})\n", + "grid.emit_template(wide)\n", + "write_chunk(wide, \"1121121\", {0: rng.uniform(0.0, 400.0, 5_000)})\n", "\n", "# Default fit=\"raise\": the 400 m span overflows the 64 m window.\n", "try:\n", @@ -402,7 +395,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": null, "id": "00ce8977", "metadata": { "execution": { @@ -412,16 +405,7 @@ "shell.execute_reply": "2026-06-25T19:05:18.356093Z" } }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "fit='degrade_resolution' -> tensor (64, 64, 128) dtype uint32\n", - "fit='collapse_bins' -> tensor (64, 64, 32) (n_bins collapsed to fit ~10 m span)\n" - ] - } - ], + "outputs": [], "source": [ "# degrade_resolution: coarsen the bins until the full span fits in 128 bins.\n", "t, _ = next(read_tensors(wide, FIELD, bottom=0.0, top=1.0, fit=\"degrade_resolution\"))\n", @@ -429,7 +413,8 @@ "\n", "# collapse_bins on a *narrow* chunk: shrink n_bins to the smallest pow2 that fits.\n", "narrow = MemoryStore()\n", - "write_chunk(narrow, FIELD, 2, {0: rng.uniform(100.0, 110.0, 4_000)})\n", + "grid.emit_template(narrow)\n", + "write_chunk(narrow, \"1121121\", {0: rng.uniform(100.0, 110.0, 4_000)})\n", "t2, _ = next(read_tensors(narrow, FIELD, bottom=0.0, top=1.0, fit=\"collapse_bins\"))\n", "print(\"fit='collapse_bins' -> tensor\", t2.shape, \"(n_bins collapsed to fit ~10 m span)\")" ] @@ -453,7 +438,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": null, "id": "8e43b0b0", "metadata": { "execution": { @@ -463,22 +448,14 @@ "shell.execute_reply": "2026-06-25T19:05:18.374485Z" } }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "chunk 42, cell 7: recovered [1. 2. 3. 4. 5.]\n", - "matches sorted input: True\n" - ] - } - ], + "outputs": [], "source": [ "from zagg.readers import read_raw_values\n", "\n", "lossless = MemoryStore()\n", + "grid.emit_template(lossless)\n", "samples = np.array([3.0, 1.0, 2.0, 5.0, 4.0]) # 5 values, well under delta -> no merges\n", - "write_chunk(lossless, FIELD, 42, {7: samples}, delta=512)\n", + "write_chunk(lossless, \"1121121\", {7: samples}, delta=512)\n", "\n", "for morton, cell_id, values in read_raw_values(lossless, FIELD):\n", " print(f\"chunk {morton}, cell {cell_id}: recovered {values}\")\n", @@ -487,7 +464,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": null, "id": "b02758c0", "metadata": { "execution": { @@ -497,19 +474,12 @@ "shell.execute_reply": "2026-06-25T19:05:18.400150Z" } }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "merged digest -> cell 0 (chunk 1) has merged centroids (weight > 1); raw values are not losslessly recoverable\n" - ] - } - ], + "outputs": [], "source": [ "# A cell with many samples at a small delta forces merges -> not recoverable.\n", "merged = MemoryStore()\n", - "write_chunk(merged, FIELD, 1, {0: rng.standard_normal(5_000)}, delta=64)\n", + "grid.emit_template(merged)\n", + "write_chunk(merged, \"1121121\", {0: rng.standard_normal(5_000)}, delta=64)\n", "try:\n", " list(read_raw_values(merged, FIELD))\n", "except ValueError as exc:\n", @@ -528,10 +498,12 @@ " flag (`uint16`/`uint32`/`float32`).\n", "- `read_raw_values(store, field)` -> exact samples per cell when the digest is\n", " merge-free, raising otherwise.\n", + "- `read_cell(store, field, cell)` -> one cell's digest by global position.\n", "\n", - "Both read the CSR ragged t-digest field zagg writes and recover the chunk morton\n", - "id from the store. Point them at a real on-disk or S3 Zarr product and the same\n", - "calls apply -- only the `store` argument changes.\n", + "All read the ragged (variable-length-bytes) t-digest field zagg writes (issue\n", + "#209) and recover the chunk morton id from the sibling `morton` coordinate.\n", + "Point them at a real on-disk or S3 Zarr product and the same calls apply -- only\n", + "the `store` argument changes.\n", "\n", "> **Note:** this notebook uses a synthetic in-memory store so it stays\n", "> Binder-runnable with no credentials. A companion check against a non-synthetic\n", @@ -560,4 +532,4 @@ }, "nbformat": 4, "nbformat_minor": 5 -} +} \ No newline at end of file diff --git a/src/zagg/config.py b/src/zagg/config.py index 66aad2f..58c38d5 100644 --- a/src/zagg/config.py +++ b/src/zagg/config.py @@ -512,6 +512,21 @@ def validate_config(config: PipelineConfig) -> None: f"accept a 'locations' keyword, which 'location' requires " f"(use a located reducer, e.g. zagg.stats.tdigest.build_tdigest)" ) + # The location channel is stored as a SIBLING vlen array named + # ``{name}_locations`` (issue #209) — a name the user never wrote. + # A field or coordinate already claiming it would silently lose in + # the template members dict and interleave into the same object + # slab at write time (data corruption), so reject at load + # (review, PR #211). + from zagg.grids.base import ragged_locations_name + + sibling = ragged_locations_name(name) + if sibling in agg_vars or sibling in config.aggregation.get("coordinates", {}): + raise ValueError( + f"Variable '{name}': its location channel is stored in a sibling " + f"array named '{sibling}', which collides with the declared " + f"field/coordinate of that name — rename one of them (issue #209)" + ) # Required per-variable keys for a temporal/event aggregation spec. ``mask`` @@ -644,7 +659,7 @@ def _validate_chunk_precompute(aggregation: dict, ds_vars: set[str]) -> None: ) from e -# Recognized per-field output kinds. ``ragged`` (CSR) is the Tier-2 carrier +# Recognized per-field output kinds. ``ragged`` is the Tier-2 carrier # for variable-length per-cell outputs; see issue #48. OUTPUT_KINDS = ("scalar", "vector", "ragged") @@ -662,7 +677,7 @@ def _validate_output_kind(name: str, meta: dict) -> None: ``ragged``). ``scalar`` fields need neither and stay the default path. ``vector`` and ``ragged`` fields may be driven by either ``function`` or ``expression``; ``len``/``count`` are rejected for both (they short-circuit - to a scalar count). See issue #29 (vector) and issue #48 (ragged/CSR). + to a scalar count). See issue #29 (vector) and issue #48 (ragged). A field may also declare ``resolution`` (``cell`` default, or ``chunk``). A ``resolution: chunk`` field (issue #30 item 2) is written ONCE per chunk @@ -672,7 +687,7 @@ def _validate_output_kind(name: str, meta: dict) -> None: ``scalar``, ``vector``, and ``ragged`` kinds may all be ``resolution: chunk`` (issue #82): a ``scalar`` companion is a plain chunk-grid array, a ``vector`` companion appends the field's ``trailing_shape`` to the chunk grid (chunked - whole), and a ``ragged`` companion is CSR at chunk granularity — one + whole), and a ``ragged`` companion holds one variable-length payload per chunk, written by ``write_ragged_to_zarr`` (phase 4c). The shape keys are validated below exactly as for cell resolution — the chunk axis just replaces the cell axis. @@ -697,7 +712,7 @@ def _validate_output_kind(name: str, meta: dict) -> None: ) # ``location`` (issue #87) is the ragged location channel; other kinds have - # no companion CSR array to carry it. + # no companion ragged array to carry it. if "location" in meta and kind != "ragged": raise ValueError( f"Variable '{name}': 'location' is only valid for kind 'ragged', not '{kind}'" @@ -715,9 +730,9 @@ def _validate_output_kind(name: str, meta: dict) -> None: raise ValueError( f"Variable '{name}': resolution '{resolution}' is not supported (allowed: {allowed})" ) - # ``ragged`` at chunk resolution (issue #82) is CSR at chunk granularity: one - # variable-length payload per chunk instead of per cell. It rides the same CSR - # writer as cell-resolution ragged (``write_ragged_to_zarr``), which collapses + # ``ragged`` at chunk resolution (issue #82) stores one variable-length + # payload per chunk instead of per cell. It rides the same vlen writer as + # cell-resolution ragged (``write_ragged_to_zarr``), which collapses # the populated cells to the single chunk payload under the same chunk-uniform # contract as scalar/vector chunk companions (raise if populated cells # disagree). No special shape key is needed beyond ``inner_shape``. @@ -774,7 +789,7 @@ def _validate_output_kind(name: str, meta: dict) -> None: # Optional location channel (issue #87): the reducer also receives the named # per-observation morton column (``locations=`` kwarg) and returns a - # ``(payload, locations)`` pair, stored as a uint64 companion CSR array. Only + # ``(payload, locations)`` pair, stored as a uint64 companion vlen array. Only # a ``function`` reducer can accept the kwarg, and the chunk-uniform collapse # of ``resolution: chunk`` has no location fold — reject both combinations. location = meta.get("location") @@ -1353,7 +1368,7 @@ def output_field_signature(config: PipelineConfig) -> list[dict]: "inner_shape": list(sig["inner_shape"]), "dtype": sig["dtype"], } - # A location channel changes the store schema (a uint64 companion CSR + # A location channel changes the store schema (a uint64 companion vlen # array — issue #87), so it belongs in the signature; keyed only when set # so existing shard-map signatures are byte-identical. if sig["location"] is not None: diff --git a/src/zagg/configs/atl03_gain_bias_healpix.yaml b/src/zagg/configs/atl03_gain_bias_healpix.yaml index ba9366e..88acdfc 100644 --- a/src/zagg/configs/atl03_gain_bias_healpix.yaml +++ b/src/zagg/configs/atl03_gain_bias_healpix.yaml @@ -124,7 +124,7 @@ output: chunk_inner: 13 # inner Zarr chunk: 4^(19-13) = 64x64 cells (K=16 per shard) # ShardingCodec output (issue #108/#133, decided on the #170 read-matrix # evidence 2026-07-07): at identical reads, sharded collapses the - # per-inner-chunk aggregation grouping + K CSR writes (o9: agg 165->19 s, + # per-inner-chunk aggregation grouping + K per-inner-chunk writes (o9: agg 165->19 s, # write 148->15 s per worker; -44% cost vs inner). The benchmark matrix # keeps both codecs as explicit A/B columns. sharded: true diff --git a/src/zagg/configs/atl03_tdigest_healpix.yaml b/src/zagg/configs/atl03_tdigest_healpix.yaml index ce657b7..872eac6 100644 --- a/src/zagg/configs/atl03_tdigest_healpix.yaml +++ b/src/zagg/configs/atl03_tdigest_healpix.yaml @@ -1,7 +1,7 @@ # ATL03 photon-height t-digest template (HEALPix, multi-chunk-per-worker). # # Per-cell approximate-quantile sketch: each grid cell stores a t-digest of its -# photon heights (h_ph) as a ``kind: ragged`` CSR field. A t-digest is a small, +# photon heights (h_ph) as a ``kind: ragged`` vlen-bytes field. A t-digest is a small, # mergeable set of (mean, weight) centroids from which any quantile (median, # p5/p95, IQR, ...) is recovered at read time -- far cheaper than storing every # photon, and exact-enough in the tails. The reducer is the pure-numpy @@ -14,9 +14,11 @@ # resolves past the old order-18 reference cap). Multi-chunk-per-worker (issue # #30 item 3) is on: one shard (the dispatch unit, parent_order 11 = 256x256 # cells) owns K = 16 finer Zarr chunks (chunk_inner 13 = 64x64 cells each). The -# worker reads the shard's granules once and writes 16 chunk regions + 16 CSR -# groups, one per inner chunk. chunk_inner unset would collapse to K=1 (shard == -# chunk), byte-identical to the single-chunk path. +# worker reads the shard's granules once and writes 16 chunk regions plus the +# ragged t-digest as one vlen-bytes object per shard (issue #209 -- the sharded +# vlen-bytes layout replacing the old per-inner-chunk CSR groups). chunk_inner +# unset would collapse to K=1 (shard == chunk), byte-identical to the +# single-chunk path. # # Confidence filter, multi-level form and read_plan match atl03.yaml (the planned # segment->photon IO of issue #43 applies here too); only the aggregation differs. @@ -67,11 +69,11 @@ aggregation: dtype: int32 fill_value: 0 h_tdigest: - # Per-cell t-digest of photon heights: a ragged (CSR) field whose payload + # Per-cell t-digest of photon heights: a ragged (vlen-bytes) field whose payload # is an (n_centroids, 2) array of (mean, weight) centroids. ``build_tdigest`` # is called as build_tdigest(h_ph_values, delta=256); ``inner_shape: [2]`` # declares the per-element centroid width. To store ONE digest per chunk - # instead of per cell, add ``resolution: chunk`` (one CSR payload per inner + # instead of per cell, add ``resolution: chunk`` (one payload per inner # chunk under the chunk-uniform contract). kind: ragged function: "zagg.stats.tdigest.build_tdigest" @@ -90,7 +92,7 @@ output: chunk_inner: 13 # inner Zarr chunk: 4^(19-13) = 64x64 cells (K=16 per shard) # ShardingCodec output (issue #108/#133, decided on the #170 read-matrix # evidence 2026-07-07): at identical reads, sharded collapses the - # per-inner-chunk aggregation grouping + K CSR writes (o9: agg 165->19 s, + # per-inner-chunk aggregation grouping + K per-inner-chunk writes (o9: agg 165->19 s, # write 148->15 s per worker; -44% cost vs inner). The benchmark matrix # keeps both codecs as explicit A/B columns. sharded: true diff --git a/src/zagg/configs/atl03_tdigest_located_healpix.yaml b/src/zagg/configs/atl03_tdigest_located_healpix.yaml index 0ef9077..7a872ce 100644 --- a/src/zagg/configs/atl03_tdigest_located_healpix.yaml +++ b/src/zagg/configs/atl03_tdigest_located_healpix.yaml @@ -6,8 +6,8 @@ # HEALPix cell enclosing the photons merged into it, computed by folding their # order-29 point-kind morton words (``HealpixGrid.assign``) through # ``mortie.common_ancestor``. A one-photon centroid stores its exact order-29 -# point word. On disk the locations are a uint64 companion CSR array -# (``{field}/{shard}/locations``) sharing offsets/cell_ids with the float32 +# point word. On disk the locations are a uint64 companion vlen array +# (the ``{field}_locations`` sibling, issue #209) row-aligned with the float32 # digest -- one logical multivariable field, two physical arrays, each channel # lossless in its own dtype (an order-29 word exceeds float64's 2^53 exact-int # ceiling, so no interleaved float layout can carry it). Read the channel back @@ -64,7 +64,7 @@ aggregation: # per-photon order-29 morton words (``locations=leaf_id``) and returns a # ``(digest, locations)`` pair -- the (n_centroids, 2) float32 (mean, # weight) array plus the (n_centroids,) uint64 per-centroid location - # vector, stored as the companion CSR array described above. + # vector, stored as the companion vlen array described above. kind: ragged function: "zagg.stats.tdigest.build_tdigest" source: h_ph @@ -83,7 +83,7 @@ output: chunk_inner: 13 # inner Zarr chunk: 4^(19-13) = 64x64 cells (K=16 per shard) # ShardingCodec output (issue #108/#133, decided on the #170 read-matrix # evidence 2026-07-07): at identical reads, sharded collapses the - # per-inner-chunk aggregation grouping + K CSR writes (o9: agg 165->19 s, + # per-inner-chunk aggregation grouping + K per-inner-chunk writes (o9: agg 165->19 s, # write 148->15 s per worker; -44% cost vs inner). The benchmark matrix # keeps both codecs as explicit A/B columns. sharded: true diff --git a/src/zagg/csr.py b/src/zagg/csr.py deleted file mode 100644 index 14493e4..0000000 --- a/src/zagg/csr.py +++ /dev/null @@ -1,266 +0,0 @@ -"""CSR (Compressed Sparse Row) writer and reader for ragged Zarr v3 arrays. - -A ragged field (issue #48) stores variable-length per-cell payloads using three -co-located Zarr v3 arrays under a common prefix:: - - {field_name}/values -- flat concatenation of all per-cell arrays - {field_name}/offsets -- per-populated-cell start index into values (length - n_populated + 1; offsets[-1] == len(values)) - {field_name}/cell_ids -- which cells (by position in the chunk's children - array) have data (length n_populated) - -This mirrors the standard CSR layout: ``values[offsets[k]:offsets[k+1]]`` is the -payload for cell ``cell_ids[k]``. Cells absent from ``cell_ids`` have no data. - -A *located* ragged field (issue #87) adds a fourth array:: - - {field_name}/locations -- flat uint64 morton location words, one per - values row, SHARING offsets/cell_ids - -so ``locations[offsets[k]:offsets[k+1]]`` is cell ``cell_ids[k]``'s per-element -location vector. The shared offsets are valid because :func:`write_csr` -enforces equal per-cell lengths. Value-only fields are byte-identical to the -pre-#87 three-array layout (no probe, no extra key). - -Profiling note (IO vs compute) -------------------------------- -For a typical ATL06 shard (4096 cells, ~100–1000 t-digest centroids per cell), -the flat ``values`` array is small relative to the scalar arrays written by the -dense path, so CSR IO overhead is dominated by the additional per-shard Zarr -metadata writes (3 arrays instead of 1). The per-shard compute for the CSR -pack/unpack (``np.concatenate`` + ``np.cumsum``) is O(total_elements) and -negligible vs t-digest construction. At δ=512, each cell stores up to 2048 -floats (1024 centroids × 2 columns); at the shard scale the flat values array -is ≲8 MB, fitting comfortably in a single Zarr chunk with ``chunks=None`` -(store as one object). For production the chunk size should be tuned to the -expected number of populated cells × centroids_per_cell × item_size. -""" - -from __future__ import annotations - -import numpy as np -import zarr -from zarr.abc.store import Store - - -def write_csr( - store: Store, - field_name: str, - values_list: list[np.ndarray], - cell_ids: list[int], - *, - dtype: str | np.dtype = "float32", - zarr_format: int = 3, - locations_list: list[np.ndarray] | None = None, -) -> None: - """Write a ragged field to a Zarr store as three CSR arrays. - - Parameters - ---------- - store : Store - Zarr-compatible store (memory, local, S3, …). - field_name : str - Prefix under which the three arrays are written: - ``{field_name}/values``, ``{field_name}/offsets``, - ``{field_name}/cell_ids``. - values_list : list of ndarray - Per-populated-cell payloads. Each element is a numpy array of shape - ``(n_elements, *inner_shape)`` (or 1-D ``(n_elements,)`` for a scalar - inner type). Empty arrays are silently skipped. - cell_ids : list of int - Cell index (position in the chunk's ``children`` array) for each entry - in ``values_list``. Must have the same length as ``values_list``. - dtype : str or numpy dtype, optional - Dtype for the ``values`` array. Default ``"float32"``. - zarr_format : int, optional - Zarr format version. Default 3. - locations_list : list of ndarray, optional - Per-cell ``uint64`` morton location vectors (issue #87), index-aligned - with ``values_list``; each element has one word per payload element. - Written as a fourth array ``{field_name}/locations`` **sharing** - ``offsets``/``cell_ids`` with ``values`` (the shared offsets are valid - because the per-cell lengths are enforced equal). ``None`` (default) - writes exactly the pre-#87 three arrays, byte-identical. - - Raises - ------ - ValueError - If ``values_list`` and ``cell_ids`` have different lengths, if any - element of ``values_list`` has an inconsistent inner shape, or if a - ``locations_list`` entry's length disagrees with its payload. - """ - if len(values_list) != len(cell_ids): - raise ValueError( - f"values_list (len {len(values_list)}) and cell_ids (len {len(cell_ids)}) " - "must have the same length" - ) - if locations_list is not None and len(locations_list) != len(values_list): - raise ValueError( - f"locations_list (len {len(locations_list)}) and values_list " - f"(len {len(values_list)}) must have the same length" - ) - - dtype = np.dtype(dtype) - - # Filter out empty payloads (cells with no data contribute nothing to CSR). - # Locations ride the same filter so they stay index-aligned. Check the empty - # entries' locations BEFORE filtering — a non-empty location vector on an - # empty payload is an upstream misalignment and must raise, not vanish. - locs_or_none: list[np.ndarray | None] = ( - list(locations_list) if locations_list is not None else [None] * len(values_list) - ) - if locations_list is not None: - for i, (arr, loc) in enumerate(zip(values_list, locations_list)): - if np.asarray(arr).size == 0 and np.asarray(loc).size > 0: - raise ValueError( - f"locations at index {i} are non-empty but the payload is empty; " - f"the location channel must stay aligned with the payload" - ) - non_empty = [ - (arr, cid, loc) - for arr, cid, loc in zip(values_list, cell_ids, locs_or_none) - if np.asarray(arr).size > 0 - ] - - flat_locs: np.ndarray | None = None - if not non_empty: - flat = np.empty(0, dtype=dtype) - offsets = np.zeros(1, dtype=np.int64) - ids_arr = np.empty(0, dtype=np.int64) - if locations_list is not None: - flat_locs = np.empty(0, dtype=np.uint64) - else: - arrays, ids, locs = zip(*non_empty) - arrays = [np.asarray(a, dtype=dtype) for a in arrays] - # Validate consistent inner shape. - inner = arrays[0].shape[1:] if arrays[0].ndim > 1 else () - for i, a in enumerate(arrays): - a_inner = a.shape[1:] if a.ndim > 1 else () - if a_inner != inner: - raise ValueError( - f"Inconsistent inner shape at index {i}: expected {inner}, got {a_inner}" - ) - flat = np.concatenate([a.reshape(a.shape[0], -1) if a.ndim > 1 else a for a in arrays]) - lengths = np.array([a.shape[0] for a in arrays], dtype=np.int64) - offsets = np.concatenate([[0], np.cumsum(lengths)]) - ids_arr = np.array(ids, dtype=np.int64) - if locations_list is not None: - locs = [np.asarray(loc, dtype=np.uint64) for loc in locs] - for i, (a, loc) in enumerate(zip(arrays, locs)): - if loc.shape != (a.shape[0],): - raise ValueError( - f"locations at index {i} have shape {loc.shape}, expected " - f"({a.shape[0]},) to share the payload's offsets" - ) - flat_locs = np.concatenate(locs) if locs else np.empty(0, dtype=np.uint64) - - _write_array(store, f"{field_name}/values", flat, zarr_format=zarr_format) - _write_array(store, f"{field_name}/offsets", offsets, zarr_format=zarr_format) - _write_array(store, f"{field_name}/cell_ids", ids_arr, zarr_format=zarr_format) - if flat_locs is not None: - _write_array(store, f"{field_name}/locations", flat_locs, zarr_format=zarr_format) - - -def _write_array(store: Store, path: str, data: np.ndarray, *, zarr_format: int = 3) -> None: - """Write a 1-D or 2-D numpy array to a Zarr store at ``path``.""" - arr = zarr.open_array( - store, - path=path, - mode="w", - shape=data.shape, - chunks=data.shape - if data.size > 0 - else tuple(max(1, d) for d in data.shape) - if data.ndim > 1 - else (1,), - dtype=data.dtype, - zarr_format=zarr_format, - ) - if data.size > 0: - arr[...] = data - - -def read_csr( - store: Store, - field_name: str, - *, - zarr_format: int = 3, - locations: bool = False, -) -> dict[str, np.ndarray]: - """Read CSR arrays written by :func:`write_csr`. - - Parameters - ---------- - store : Store - Zarr-compatible store. - field_name : str - Prefix used when writing (``{field_name}/values`` etc.). - zarr_format : int, optional - Zarr format version. Must match the version used when writing. - Default 3. - locations : bool, optional - Also read the ``{field_name}/locations`` companion (issue #87) into a - ``"locations"`` key. Opt-in so value-only reads stay free of the extra - store probe. Raises a clear error when the field carries no location - channel. - - Returns - ------- - dict with keys ``"values"``, ``"offsets"``, ``"cell_ids"``. - - ``values`` : flat concatenation of all per-cell payloads. - - ``offsets`` : length ``n_populated + 1``; ``values[offsets[k]:offsets[k+1]]`` - is the payload for cell ``cell_ids[k]``. - - ``cell_ids`` : which cells (by chunk position) have data. - With ``locations=True``, also ``"locations"`` — flat uint64 morton words - sharing ``offsets`` with ``values``. - """ - values = np.asarray( - zarr.open_array(store, path=f"{field_name}/values", mode="r", zarr_format=zarr_format)[...] - ) - offsets = np.asarray( - zarr.open_array(store, path=f"{field_name}/offsets", mode="r", zarr_format=zarr_format)[...] - ) - cell_ids = np.asarray( - zarr.open_array(store, path=f"{field_name}/cell_ids", mode="r", zarr_format=zarr_format)[ - ... - ] - ) - out = {"values": values, "offsets": offsets, "cell_ids": cell_ids} - if locations: - try: - out["locations"] = np.asarray( - zarr.open_array( - store, path=f"{field_name}/locations", mode="r", zarr_format=zarr_format - )[...] - ) - except (FileNotFoundError, KeyError) as e: - raise ValueError( - f"{field_name!r} has no locations array; it was not written as a " - f"located ragged field (declare 'location:' on the field — issue #87)" - ) from e - return out - - -def iter_csr_cells( - csr: dict[str, np.ndarray], -) -> list[tuple[int, np.ndarray]]: - """Decode a CSR dict into a list of ``(cell_id, payload)`` pairs. - - Parameters - ---------- - csr : dict - As returned by :func:`read_csr`. - - Returns - ------- - list of (int, ndarray) - One entry per populated cell, in the order they were written. - """ - values = csr["values"] - offsets = csr["offsets"] - cell_ids = csr["cell_ids"] - result = [] - for k, cid in enumerate(cell_ids): - payload = values[offsets[k] : offsets[k + 1]] - result.append((int(cid), payload)) - return result diff --git a/src/zagg/grids/base.py b/src/zagg/grids/base.py index a0d9a1c..ff06090 100644 --- a/src/zagg/grids/base.py +++ b/src/zagg/grids/base.py @@ -23,6 +23,8 @@ from __future__ import annotations +import warnings +from contextlib import contextmanager from typing import Any, Protocol, runtime_checkable import numpy as np @@ -30,6 +32,167 @@ ShardKey = Any # int for HEALPix, tuple[int,int] for rectilinear, etc. +#: Array-attrs key on a ragged vlen-bytes array (issue #209) recording the +#: per-cell element interpretation. The value is +#: ``{"element": {"dtype": "", "shape": [-1, *inner_shape]}}``: +#: each populated cell's value is the raw little-endian bytes of an +#: ``(n, *inner_shape)`` array (``-1`` marks the per-cell varying count), so a +#: reader reconstructs cell ``i`` as +#: ``np.frombuffer(a[i], dtype).reshape(-1, *inner_shape)``. A LOCATED field's +#: payload array additionally carries ``{"locations": ""}`` +#: (issue #87) — the reader binds the uint64 channel by that declaration, not +#: by reconstructing the naming convention (review, PR #211). The block is +#: versioned (``{"spec": RAGGED_SPEC}``, the coverage-envelope discipline) so +#: readers fail loudly on a future revision instead of half-parsing it. +RAGGED_ELEMENT_ATTR = "ragged" + +#: Convention version stamped into the :data:`RAGGED_ELEMENT_ATTR` block and +#: strict-checked by the readers (``readers/tdigest_tensor._open_ragged``). +#: This attrs seam is the INTERIM contract: the issue #210 typed +#: ``vlen-array`` dtype migration moves the element declaration into the +#: zarr data type itself and supersedes (bumps or removes) this marker. +RAGGED_SPEC = "zagg-ragged/1" + +#: zstd level of the ragged inner codec chain — 3, matching the coverage +#: sidecar precedent (``zagg.hive._ZSTD_LEVEL``), fixed so identical payloads +#: produce identical objects across workers. +RAGGED_ZSTD_LEVEL = 3 + + +def ragged_locations_name(field_name: str) -> str: + """On-disk array name of a located ragged field's uint64 channel (issue #87). + + Under the vlen-bytes layout (issue #209) the location words are a SIBLING + vlen array (``{field}_locations``) row-aligned with the digest payload — + the CSR layout's fourth in-group array cannot nest under what is now an + array node. + """ + return f"{field_name}_locations" + + +@contextmanager +def vlen_dtype_warning_suppressed(): + """Suppress zarr's vlen-bytes dtype-naming warning at array CREATION only. + + zarr-python names the dtype ``variable_length_bytes`` in metadata while + the v3 registry name is ``bytes`` (zarr-python#3517, accepted both ways on + read), so creating a ragged vlen array emits an + ``UnstableSpecificationWarning`` about that naming. Scoped to the exact + message (the ``coverage.moc`` suppression precedent) so nothing else is + silenced; reads/writes/opens do not warn. + """ + from zarr.errors import UnstableSpecificationWarning + + with warnings.catch_warnings(): + warnings.filterwarnings( + "ignore", + message=r"The data type \(VariableLengthBytes\(\)\)", + category=UnstableSpecificationWarning, + ) + yield + + +def ragged_array_spec( + *, + shape, + dims, + inner_chunk_shape, + shard_shape=None, + element_dtype, + inner_shape=(), + locations=None, +): + """Vlen-bytes ``ArraySpec`` for a ``kind: ragged`` field (issue #209). + + The sharded vlen-bytes layout: ONE zarr array with the + ``variable_length_bytes`` data type replaces the per-inner-chunk CSR + subgroups (~7 objects per populated inner chunk). Each populated cell + holds the raw little-endian bytes of its ``(n, *inner_shape)`` payload + (``n`` varies per cell); empty cells keep the ``b""`` fill, and an + all-empty inner chunk is omitted from the shard index — the same + sub-shard sparsity the dense arrays get. The element interpretation is + self-describing via :data:`RAGGED_ELEMENT_ATTR` in the array attrs. + + Codec chain: ``[vlen-bytes, zstd(level=3)]``. The zstd deviates from the + dense arrays' bytes-only/uncompressed policy deliberately: a vlen payload + has no fixed-width raw layout to preserve, and level 3 matches the + coverage-sidecar precedent. With ``shard_shape`` the chain rides INSIDE a + ``ShardingCodec`` (outer chunk == ``shard_shape``), collapsing a shard's K + inner chunks to one object with an internal index — single-cell reads stay + 2 GETs (index suffix + one ranged inner chunk). ``None`` keeps a regular + array chunked at ``inner_chunk_shape`` (one object per inner chunk — the + unsharded per-chunk-write layout). + + Parameters + ---------- + shape : tuple of int + Array shape — the grid's cell axes (or the chunk grid for a + ``resolution: chunk`` companion). + dims : tuple of str + Dimension names matching ``shape``. + inner_chunk_shape : tuple of int + Read-chunk shape (``grid.chunk_shape``; ``(1,)*ndim`` for a chunk + companion). + shard_shape : tuple of int, optional + ShardingCodec outer chunk. ``None`` (default) emits a regular array. + element_dtype : str + Numpy dtype of one payload element (recorded in attrs; the bytes are + little-endian). + inner_shape : tuple of int, optional + Per-element trailing shape (``sig["inner_shape"]``, e.g. ``(2,)`` for + a centroid pair). Empty for a flat per-cell vector. + locations : str, optional + Name of the located field's uint64 sibling array (issue #87), + declared in the payload array's attrs so a reader binds the channel + by METADATA, not by reconstructing the naming convention (review, + PR #211). ``None`` (unlocated) records nothing. + + Returns + ------- + ArraySpec + """ + from pydantic_zarr.experimental.v3 import ArraySpec, NamedConfig + + inner_codecs = [ + {"name": "vlen-bytes", "configuration": {}}, + {"name": "zstd", "configuration": {"level": RAGGED_ZSTD_LEVEL, "checksum": False}}, + ] + if shard_shape is not None: + # Index codecs mirror ``sharded_array_spec`` (zarr's create_array default). + codecs: tuple = ( + NamedConfig( + name="sharding_indexed", + configuration={ + "chunk_shape": [int(c) for c in inner_chunk_shape], + "codecs": inner_codecs, + "index_codecs": [ + {"name": "bytes", "configuration": {"endian": "little"}}, + {"name": "crc32c"}, + ], + "index_location": "end", + }, + ), + ) + chunk_shape = tuple(int(c) for c in shard_shape) + else: + codecs = tuple(NamedConfig(**c) for c in inner_codecs) + chunk_shape = tuple(int(c) for c in inner_chunk_shape) + element = {"dtype": str(element_dtype), "shape": [-1, *(int(s) for s in inner_shape)]} + ragged_meta: dict = {"spec": RAGGED_SPEC, "element": element} + if locations is not None: + ragged_meta["locations"] = str(locations) + return ArraySpec( + attributes={RAGGED_ELEMENT_ATTR: ragged_meta}, + shape=tuple(int(s) for s in shape), + dimension_names=tuple(dims), + data_type="variable_length_bytes", + chunk_grid=NamedConfig(name="regular", configuration={"chunk_shape": list(chunk_shape)}), + chunk_key_encoding=NamedConfig(name="default", configuration={"separator": "/"}), + codecs=codecs, + storage_transformers=(), + fill_value="", + ) + def vector_array_spec(base, sig, *, base_dims, base_chunk_shape): """Extend a scalar data-var ``ArraySpec`` with a trailing payload dim. @@ -261,8 +424,8 @@ def block_index(self, shard_key: ShardKey) -> tuple[int, ...]: def shard_label(self, shard_key: ShardKey) -> str: """External string form of a shard key (issue #199). - Used wherever a shard id surfaces outside the process — CSR subgroup - names, async ``.status`` object keys, log lines. HEALPix renders the + Used wherever a shard id surfaces outside the process — hive leaf + ids, async ``.status`` object keys, log lines. HEALPix renders the packed word as its decimal morton string (D1 in ``docs/design/sparse_coverage.md``); rectilinear keeps the packed tile int's decimal digits. @@ -327,10 +490,16 @@ def shard_label(grid, shard_key) -> str: __all__ = [ "OutputGrid", + "RAGGED_ELEMENT_ATTR", + "RAGGED_SPEC", + "RAGGED_ZSTD_LEVEL", "ShardKey", "InconsistentShardError", "shard_label", + "ragged_array_spec", + "ragged_locations_name", "vector_array_spec", + "vlen_dtype_warning_suppressed", "chunk_array_spec", "sharded_array_spec", ] diff --git a/src/zagg/grids/healpix.py b/src/zagg/grids/healpix.py index 427065d..ef678d1 100644 --- a/src/zagg/grids/healpix.py +++ b/src/zagg/grids/healpix.py @@ -21,8 +21,11 @@ from zagg.grids.base import ( InconsistentShardError, chunk_array_spec, + ragged_array_spec, + ragged_locations_name, sharded_array_spec, vector_array_spec, + vlen_dtype_warning_suppressed, ) from zagg.grids.morton import morton_decimal, to_morton_array @@ -462,7 +465,7 @@ def shard_label(self, shard_key) -> str: """Decimal morton string for this shard's packed word (issue #199). The external form of a HEALPix shard id (D1 in - ``docs/design/sparse_coverage.md``): CSR subgroup names, ``.status`` + ``docs/design/sparse_coverage.md``): hive leaf ids, ``.status`` object keys, and log lines all carry e.g. ``-31123``, never the raw packed-word integer. """ @@ -573,7 +576,9 @@ def nests_with(self, other) -> bool: def emit_template(self, store: Store, *, overwrite: bool = False) -> Store: """Write the Zarr template (group + arrays) to ``store``.""" spec = self._spec() - with zarr_config.set({"async.concurrency": 128}): + # Ragged vlen-array creation warns about the dtype NAME only + # (zarr-python#3517); message-scoped suppression, see grids.base. + with zarr_config.set({"async.concurrency": 128}), vlen_dtype_warning_suppressed(): spec.to_zarr(store, self.group_path, overwrite=overwrite) return store @@ -593,7 +598,9 @@ def emit_shard_template(self, store: Store, *, overwrite: bool = False) -> Store # the single shared store. raise ValueError("hive leaf templates do not support sharded output") spec = GroupSpec(members={self.group_path: self.shard_spec()}, attributes={}) - with zarr_config.set({"async.concurrency": 128}): + # Ragged vlen-array creation warns about the dtype NAME only + # (zarr-python#3517); message-scoped suppression, see grids.base. + with zarr_config.set({"async.concurrency": 128}), vlen_dtype_warning_suppressed(): spec.to_zarr(store, "", overwrite=overwrite) return store @@ -606,9 +613,11 @@ def shard_spec(self) -> GroupSpec: Identical member set to :meth:`spec` — same dtypes, fills, chunking — with the cells axis sized to one shard and the ``resolution: chunk`` - companions sized to the shard's K inner chunks. + companions sized to the shard's K inner chunks. The one deliberate + difference (``leaf=True``): a ragged field's vlen array shards across + the whole leaf, so a shard's digest is ONE object (issue #209). """ - return self._group_spec(self.cells_per_shard, (self.chunks_per_shard,)) + return self._group_spec(self.cells_per_shard, (self.chunks_per_shard,), leaf=True) # ── internals ──────────────────────────────────────────────────────── @@ -619,7 +628,9 @@ def _spec(self) -> GroupSpec: n_pixels = self.n_children * self.n_shards return self._group_spec(n_pixels, self.chunk_grid_shape) - def _group_spec(self, n_pixels: int, chunk_grid_shape: tuple[int, ...]) -> GroupSpec: + def _group_spec( + self, n_pixels: int, chunk_grid_shape: tuple[int, ...], *, leaf: bool = False + ) -> GroupSpec: base = ArraySpec( attributes={}, shape=(n_pixels,), @@ -670,13 +681,49 @@ def _shard(arr): members["aoi_mask"] = _shard(base.with_data_type("bool").with_fill_value(False)) for name, meta in get_agg_fields(self.config).items(): sig = get_output_signature(meta) - # Ragged fields (issue #48) are stored as CSR subgroups - # (``{name}/{shard_key}/values|offsets|cell_ids``) written fresh by - # ``write_ragged_to_zarr`` — NOT a dense array. Emitting a dense - # ``{name}`` array here would make ``{name}`` an array, and CSR's - # per-shard child groups under it would then collide ("only groups may - # have child nodes"). Skip them so ``{name}`` stays a group prefix. + # Ragged fields (issue #48) are ONE vlen-bytes array on the cell + # grid (issue #209 — the sharded vlen-bytes layout replacing the + # per-inner-chunk CSR subgroups). At ``resolution: chunk`` the + # array sits on the chunk grid instead (one payload per chunk, the + # ragged analogue of the scalar/vector companions). A located + # field (issue #87) adds a sibling uint64 vlen array. if sig["kind"] == "ragged": + if sig["resolution"] == "chunk": + rag_kw: dict = { + "shape": chunk_grid_shape, + "dims": ("chunks",), + "inner_chunk_shape": (1,), + } + else: + # The ShardingCodec outer chunk mirrors the dense arrays + # when ``sharded`` (the issue #133 object split included); + # a hive LEAF (unsharded, K inner chunks in one + # self-describing store) shards the ragged field across + # the whole leaf so a shard's digest is ONE object. The + # unsharded FLAT layout stays a regular array (one object + # per inner chunk — the streaming per-chunk write must not + # read-modify-write a shared shard object). + if self.sharded: + rag_shard: tuple | None = (self.cells_per_shard_object,) + elif leaf and self.chunks_per_shard > 1: + rag_shard = (self.cells_per_shard,) + else: + rag_shard = None + rag_kw = { + "shape": (n_pixels,), + "dims": ("cells",), + "inner_chunk_shape": (self.cells_per_chunk,), + "shard_shape": rag_shard, + } + located = ragged_locations_name(name) if sig.get("location") else None + members[name] = ragged_array_spec( + element_dtype=sig["dtype"] or "float32", + inner_shape=sig["inner_shape"], + locations=located, + **rag_kw, + ) + if located: + members[located] = ragged_array_spec(element_dtype="uint64", **rag_kw) continue dtype = meta.get("dtype", "float32") fill = meta.get("fill_value", "NaN") diff --git a/src/zagg/grids/morton.py b/src/zagg/grids/morton.py index 56c389c..20941f3 100644 --- a/src/zagg/grids/morton.py +++ b/src/zagg/grids/morton.py @@ -73,7 +73,7 @@ def morton_decimal(word) -> str: The external/path form of a shard id per the sparse-coverage design record (``docs/design/sparse_coverage.md`` D1): the packed ``uint64`` word stays the canonical in-memory/wire form, and every externally visible string — - CSR subgroup names, ``.status`` object keys, log lines — renders through + hive leaf ids, ``.status`` object keys, log lines — renders through mortie's decode-through-kernel decimal repr (e.g. ``-31123``). Raises ``ValueError`` on an empty, invalid, or negative word (a path component must never be silently wrong) — a NEGATIVE int here is usually a *legacy @@ -98,7 +98,7 @@ def morton_word(label: str) -> int: """Parse a decimal morton string back to its packed word (issue #199). The inverse of :func:`morton_decimal` at the zagg boundary — used where an - external decimal id re-enters (``--morton-cell``, CSR subgroup names on the + external decimal id re-enters (``--morton-cell``, hive leaf ids on the read path). Raises ``ValueError`` on a malformed id. Implementation note: this rides mortie's private-but-documented diff --git a/src/zagg/grids/rectilinear.py b/src/zagg/grids/rectilinear.py index 89be176..351e43e 100644 --- a/src/zagg/grids/rectilinear.py +++ b/src/zagg/grids/rectilinear.py @@ -50,7 +50,14 @@ get_output_signature, output_field_signature, ) -from zagg.grids.base import chunk_array_spec, sharded_array_spec, vector_array_spec +from zagg.grids.base import ( + chunk_array_spec, + ragged_array_spec, + ragged_locations_name, + sharded_array_spec, + vector_array_spec, + vlen_dtype_warning_suppressed, +) OOB_SENTINEL: int = -1 @@ -520,7 +527,9 @@ def emit_template(self, store: Store, *, overwrite: bool = False) -> Store: from zarr import open_array spec = self._spec() - with zarr_config.set({"async.concurrency": 128}): + # Ragged vlen-array creation warns about the dtype NAME only + # (zarr-python#3517); message-scoped suppression, see grids.base. + with zarr_config.set({"async.concurrency": 128}), vlen_dtype_warning_suppressed(): spec.to_zarr(store, self.group_path, overwrite=overwrite) # Populate the x/y coord arrays with cell-centre coordinates so # downstream readers (xarray, rioxarray) get usable spatial axes. @@ -596,11 +605,34 @@ def _shard(arr): members["aoi_mask"] = _shard(base.with_data_type("bool").with_fill_value(False)) for name, meta in get_agg_fields(self.config).items(): sig = get_output_signature(meta) - # Ragged fields (issue #48) are CSR subgroups written fresh by - # ``write_ragged_to_zarr`` (``{name}/{shard_key}/...``), not a dense - # array — skip them so ``{name}`` stays a group prefix and the CSR - # child nodes don't collide with a dense array at the same path. + # Ragged fields (issue #48) are ONE vlen-bytes array on the (y, x) + # cell grid (issue #209 — the sharded vlen-bytes layout replacing + # the per-inner-chunk CSR subgroups); ``resolution: chunk`` sits on + # the chunk grid instead. Mirrors the HEALPix branch; a located + # field (issue #87) adds a sibling uint64 vlen array. if sig["kind"] == "ragged": + if sig["resolution"] == "chunk": + rag_kw: dict = { + "shape": self.chunk_grid_shape, + "dims": ("chunk_y", "chunk_x"), + "inner_chunk_shape": (1, 1), + } + else: + rag_kw = { + "shape": self.array_shape, + "dims": ("y", "x"), + "inner_chunk_shape": self.chunk_shape, + "shard_shape": (self.chunk_h, self.chunk_w) if self.sharded else None, + } + located = ragged_locations_name(name) if sig.get("location") else None + members[name] = ragged_array_spec( + element_dtype=sig["dtype"] or "float32", + inner_shape=sig["inner_shape"], + locations=located, + **rag_kw, + ) + if located: + members[located] = ragged_array_spec(element_dtype="uint64", **rag_kw) continue dtype = meta.get("dtype", "float32") fill = meta.get("fill_value", "NaN") diff --git a/src/zagg/hive.py b/src/zagg/hive.py index 855d0e9..a76f8a5 100644 --- a/src/zagg/hive.py +++ b/src/zagg/hive.py @@ -685,21 +685,25 @@ def process_and_write_hive( The SHARED per-shard write path for both backends (phase 3): the local runner's ``_cell_work`` and the Lambda handler's hive branch both call - this, so leaf templating, chunk placement, CSR naming, and stamp ordering - cannot drift between dispatchers. The shard's output is a self-describing - leaf zarr at :func:`shard_leaf_path` ``(store_root, shard_key)`` (D3), with - dense chunks written at leaf-LOCAL block indices and — as the shard's - FINAL write — the D4 commit stamp on the leaf's root group. The leaf - template is emitted lazily on the first chunk write (mirroring the Lambda - handler's lazy store open), so a no-data shard never creates the - ``.zarr/`` prefix; a worker that dies mid-shard leaves an UNSTAMPED prefix - — debris, overwritten wholesale on retry (``overwrite=True`` on the leaf - template makes the retry idempotent). ``profile`` forwards to - ``process_shard`` (issue #100); the write phase is interleaved with the - stream on this path, so no separate ``write`` timing is recorded. + this, so leaf templating, chunk placement, ragged layout, and stamp + ordering cannot drift between dispatchers. The shard's output is a + self-describing leaf zarr at :func:`shard_leaf_path` ``(store_root, + shard_key)`` (D3), with dense chunks written at leaf-LOCAL block indices + and — as the shard's FINAL write — the D4 commit stamp on the leaf's root + group. The leaf template is emitted lazily on the first chunk write + (mirroring the Lambda handler's lazy store open), so a no-data shard never + creates the ``.zarr/`` prefix; a worker that dies mid-shard leaves an + UNSTAMPED prefix — debris, overwritten wholesale on retry + (``overwrite=True`` on the leaf template makes the retry idempotent). + ``profile`` forwards to ``process_shard`` (issue #100); the write phase is + interleaved with the stream on this path, so no separate ``write`` timing + is recorded. """ - from zagg.grids.base import shard_label - from zagg.processing import process_shard, write_dataframe_to_zarr, write_ragged_to_zarr + from zagg.processing import ( + process_shard, + write_dataframe_to_zarr, + write_ragged_leaf_to_zarr, + ) from zagg.store import open_store leaf_path = shard_leaf_path(store_root, shard_key) @@ -721,17 +725,30 @@ def _leaf(): box["store"] = store return box["store"] - single_chunk = int(getattr(grid, "chunks_per_shard", 1)) == 1 + # Ragged fields accumulate across the streamed chunks (leaf-LOCAL blocks) + # and are written ONCE after the stream (issue #209): the leaf's ragged + # vlen array is a single ShardingCodec object spanning the shard, so a + # per-chunk write here would read-modify-write that object K times. + # Memory bound (review, PR #211): this re-adds an O(shard-payload) term to + # the otherwise O(chunk) streaming path (issue #91) — at the o8 t-digest + # scale this fix exists to unlock (sparse NEON o8: 17.6 M centroids × 8 B + # ≈ 141 MB of held payload; ~200 MB peak through the single write, once + # per-cell ``bytes``-object overhead and the assembled ~60 MB shard object + # are counted). Accepted deliberately: workers run 4 GB (issue #193), the + # dense side keeps its O(chunk) stream-and-free bound, and the + # accumulation is what deletes the ~K×7-object PUT storm that was ~1/3 of + # shard wall at CONUS scale (issue #209). If 88S-scale shards or the #148 + # streaming budget ever say otherwise, the escape valve is spilling the + # ragged field back to per-inner-chunk writes against a regular-chunked + # vlen array (the unsharded flat layout) — named here, not built. + ragged_chunks: list = [] def _write_chunk(block_index, carrier, ragged): store = _leaf() local = leaf_block_index(grid, block_index, shard_key) write_dataframe_to_zarr(carrier, store, grid=grid, chunk_idx=local) - # CSR subgroup naming inside a leaf mirrors the flat layout: the shard - # label at K==1; the LOCAL chunk ordinal at K>1 (leaf arrays are - # 1-D — hive is HEALPix-only, validated at config load). - ragged_key = shard_label(grid, shard_key) if single_chunk else int(local[0]) - write_ragged_to_zarr(ragged, store, grid=grid, shard_key=ragged_key) + if ragged: + ragged_chunks.append((local, ragged)) # Occupied-cell sink (issue #200): the worker already holds the shard's # populated cell words; collect them here to derive the stamp's coverage. @@ -755,7 +772,10 @@ def _write_chunk(block_index, carrier, ragged): # box payload rides it (zero extra requests), the exact-occupancy bitmap # sidecar is PUT just before it (issue #200 phase 2), and both inherit # its debris semantics: a torn worker's coverage never becomes visible. + # The leaf write order is pinned: dense (streamed) -> ragged (one object, + # issue #209) -> coverage sidecar -> stamp. if "store" in box and not metadata.get("error"): + write_ragged_leaf_to_zarr(ragged_chunks, box["store"], grid=grid) words = np.concatenate(occupied) if occupied else None if words is not None and words.size == 0: words = None diff --git a/src/zagg/processing/__init__.py b/src/zagg/processing/__init__.py index ba949be..633b2bb 100644 --- a/src/zagg/processing/__init__.py +++ b/src/zagg/processing/__init__.py @@ -59,13 +59,13 @@ from zagg.processing.worker import process_morton_cell, process_shard from zagg.processing.write import ( _arrow_column, - _block_index_key, _build_output, _carrier_empty, _chunk_resolution_fields, _chunk_uniform_value, _iter_carrier_columns, write_dataframe_to_zarr, + write_ragged_leaf_to_zarr, write_ragged_to_zarr, write_shard_to_zarr, ) @@ -82,6 +82,7 @@ "process_morton_cell", "process_shard", "write_dataframe_to_zarr", + "write_ragged_leaf_to_zarr", "write_ragged_to_zarr", "write_shard_to_zarr", # aggregate-stage helpers @@ -113,7 +114,6 @@ "_validate_planned_config", # write-stage helpers "_arrow_column", - "_block_index_key", "_build_output", "_carrier_empty", "_chunk_resolution_fields", diff --git a/src/zagg/processing/aggregate.py b/src/zagg/processing/aggregate.py index 3c58076..53de05d 100644 --- a/src/zagg/processing/aggregate.py +++ b/src/zagg/processing/aggregate.py @@ -329,11 +329,11 @@ def _has_vector_fields(config: PipelineConfig) -> bool: def _has_ragged_fields(config: PipelineConfig) -> bool: - """Whether any aggregation field declares a ``ragged`` (CSR) output. + """Whether any aggregation field declares a ``ragged`` output. Ragged fields (issue #48) carry variable-length per-cell payloads and are - collected separately from scalar/vector fields; they are written via the CSR - writer rather than the dense Zarr path. + collected separately from scalar/vector fields; they are written via the + ragged vlen writer rather than the dense Zarr path. """ return any( get_output_signature(meta)["kind"] == "ragged" for meta in get_agg_fields(config).values() @@ -427,8 +427,8 @@ def calculate_cell_statistics( # to a Python float; a ``kind: vector`` expression is coerced through the # same ``_coerce_field_value``/``trailing_shape``/dtype path as a vector # ``function`` field (issue #29). A ``kind: ragged`` expression (issue #48) - # returns the raw result as a numpy array — the CSR writer receives it as - # a variable-length per-cell payload. + # returns the raw result as a numpy array — the ragged writer receives it + # as a variable-length per-cell payload. if expression: if sig["kind"] == "vector": out = _eval_expression_raw(expression, cell_data) @@ -479,9 +479,10 @@ def calculate_cell_statistics( # Located ragged field (issue #87): hand the reducer the named # per-observation morton column and accept a ``(payload, locations)`` - # pair — the uint64 locations ride beside the payload into the CSR - # companion array. Only the HEALPix read path supplies ``leaf_id``, so a - # missing column is a grid/config mismatch, reported clearly. + # pair — the uint64 locations ride beside the payload into the + # ``{field}_locations`` sibling vlen array. Only the HEALPix read path + # supplies ``leaf_id``, so a missing column is a grid/config mismatch, + # reported clearly. if sig["kind"] == "ragged" and sig["location"] is not None: loc_col = sig["location"] if loc_col not in cell_data: @@ -513,7 +514,7 @@ def calculate_cell_statistics( # Scalar fields stay byte-for-byte identical to the pre-#29 path; a # declared ``vector`` field coerces to its trailing_shape (issue #29); a # ``ragged`` field (issue #48) returns a variable-length numpy array that - # the CSR writer later packs into flat + offsets + cell_ids arrays. + # the ragged writer later encodes as the cell's vlen-bytes payload. if sig["kind"] == "vector": result[name] = _coerce_field_value(out, sig) elif sig["kind"] == "ragged": @@ -531,7 +532,7 @@ def _empty_cell_value(meta: dict): ``np.nan`` otherwise. A ``vector`` field (issue #29) instead gets a full ``trailing_shape`` array filled with its schema-declared sentinel (:func:`_field_sentinel`), so empty and populated cells emit the same shape. - A ``ragged`` field (issue #48) returns an empty list ``[]`` — the CSR writer + A ``ragged`` field (issue #48) returns an empty list ``[]`` — the ragged writer handles absent cells by leaving them out of ``cell_ids``. A *located* ragged field (issue #87) returns an empty ``(payload, locations)`` pair instead, keeping the pair contract uniform for direct callers. @@ -557,7 +558,7 @@ def _coerce_field_value(value, sig: dict) -> np.ndarray: """Coerce a ``vector`` field's aggregation output to its declared signature. The field's ``function`` or ``expression`` must yield exactly - ``trailing_shape`` values (issue #29 Tier-1 fixed-width vectors; ragged/CSR + ``trailing_shape`` values (issue #29 Tier-1 fixed-width vectors; ragged is Tier 2). Returns a contiguous array of the declared dtype (default ``float32``), so every cell emits an identically-shaped slab the dense writer (phase 5) can stack. @@ -577,7 +578,7 @@ def _coerce_ragged_value(value, sig: dict) -> np.ndarray: A ragged field (issue #48) emits a variable-length array of shape ``(n_elements, *inner_shape)`` per cell. This function verifies the inner dimensions match the declared ``inner_shape`` and returns a contiguous - array of the declared dtype (default ``float32``), ready for the CSR writer. + array of the declared dtype (default ``float32``), ready for the ragged writer. Parameters ---------- @@ -683,7 +684,7 @@ def _aggregate_chunk_cells( ragged_locations, cells_with_data)``: dense fields preallocated to ``(n_cells, *trailing_shape)`` and filled per cell; ragged fields collected as ``(payloads, cell_indices)`` keyed by the cell's position in ``children`` - (the chunk-local index the CSR writer expects). ``ragged_locations`` holds + (the chunk-local index the ragged writer expects). ``ragged_locations`` holds the per-cell uint64 location vectors for located fields only (issue #87) — keyed like ``ragged_payloads`` and index-aligned with them; unlocated fields have no entry. @@ -749,7 +750,7 @@ def _aggregate_chunk_cells( # Fail fast if the value's shape disagrees with the declared # signature (a located field must deliver its pair; empty cells # are exempt as they were skipped above) — a silent mismatch - # would surface much later as a CSR length error in the writer. + # would surface much later as a length error in the ragged writer. if (key in ragged_locations) != (locations is not None): raise ValueError( f"ragged field {key!r}: located/unlocated mismatch between the " diff --git a/src/zagg/processing/streaming.py b/src/zagg/processing/streaming.py index d478744..aaddc02 100644 --- a/src/zagg/processing/streaming.py +++ b/src/zagg/processing/streaming.py @@ -126,7 +126,7 @@ class StreamingAggregator: ``len`` fields by summation, tdigest fields via ``merge_tdigests``. ``chunk_outputs`` then emits the exact ``(stats_arrays, ragged_payloads, ragged_cell_indices, cells_with_data)`` shape ``_aggregate_chunk_cells`` - returns, so the worker's carrier/CSR construction is shared verbatim. + returns, so the worker's carrier/ragged construction is shared verbatim. """ def __init__(self, config: PipelineConfig, grid, handoff: str, buffer_granules: int): diff --git a/src/zagg/processing/worker.py b/src/zagg/processing/worker.py index 708c476..ef57f88 100644 --- a/src/zagg/processing/worker.py +++ b/src/zagg/processing/worker.py @@ -129,7 +129,7 @@ def process_shard( is only a no-config safety net for direct callers. pyarrow is not used on either path. ragged_out : dict, optional - Out-param sink for ``kind: ragged`` (CSR) fields (issue #48). When a dict + Out-param sink for ``kind: ragged`` fields (issue #48). When a dict is passed, it is filled in place with ``{field_name: (values_list, cell_ids)}`` — ``values_list`` the per-populated-cell payload arrays and ``cell_ids`` their position in the chunk's ``children`` block — for the @@ -156,7 +156,7 @@ def process_shard( with one ``(block_index, carrier, ragged)`` tuple per chunk — ``block_index`` the chunk's storage block (from ``grid.iter_chunks``), ``carrier`` its dense DataFrame/Table, ``ragged`` its - ``{field: (values_list, cell_ids)}`` CSR map — for the caller to write K + ``{field: (values_list, cell_ids)}`` ragged map — for the caller to write K regions + K companion slices. The returned 2-tuple's ``df_out`` is an empty carrier in that case (the real carriers live in ``chunk_results``). ``None`` (default) is the K==1 path: the single chunk's carrier is the @@ -198,7 +198,7 @@ def process_shard( (DataFrame, metadata) DataFrame in canonical chunk order; metadata dict with ``shard_key``, ``cells_with_data``, ``total_obs``, ``granule_count``, - ``files_processed``, ``duration_s``, ``error``. Ragged (CSR) fields are + ``files_processed``, ``duration_s``, ``error``. Ragged fields are delivered out-of-band via ``ragged_out`` (above), not in this tuple. At K>1 the per-chunk carriers + ragged are delivered via ``chunk_results``. """ diff --git a/src/zagg/processing/write.py b/src/zagg/processing/write.py index 95b6f3b..ea6eaef 100644 --- a/src/zagg/processing/write.py +++ b/src/zagg/processing/write.py @@ -7,37 +7,19 @@ stays acyclic. """ -import logging -from concurrent.futures import ThreadPoolExecutor, as_completed - import numpy as np import pandas as pd from zarr import config, open_array from zarr.abc.store import Store -from zagg.concurrency import fd_safe_max_workers from zagg.config import ( PipelineConfig, get_agg_fields, get_output_signature, ) -from zagg.csr import write_csr +from zagg.grids.base import ragged_locations_name from zagg.grids.morton import is_morton_array, is_morton_arrow, morton_to_arrow, morton_words -# The sharded path's K ragged (CSR) subgroup writes fan out over a bounded thread -# pool (issue #142): each inner chunk emits an independent CSR group at a disjoint -# prefix, so they are write-independent and latency-bound (~K×3 tiny objects), and -# a serial loop dominated a t-digest shard's write time. 128 mirrors the dense -# path's ``async.concurrency`` budget (issue #108). -_RAGGED_WRITE_CONCURRENCY = 128 - -# Cap on the per-failure roll call logged when sharded CSR writes fail -# (issue #186): enough to show the race signature across subgroups without -# flooding CloudWatch when many of the K writes fail at once. -_RAGGED_FAILURE_LOG_CAP = 20 - -logger = logging.getLogger(__name__) - def _arrow_column(block: np.ndarray, sig: dict): """Build the Arrow column for one agg field from its per-cell stats block. @@ -275,9 +257,12 @@ def write_shard_to_zarr( column (data vars + coords) is placed into the shard slab by the chunk's own region (``grid.shard_local_region(block_index, shard_key)``), reshaped to the inner chunk's cell shape (``grid.chunk_shape``) — grid-agnostic, so HEALPix - (1-D) and rectilinear (2-D) share one path. The ``resolution: chunk`` - companions and ragged (CSR) fields are NOT sharded — they are written per inner - chunk via the existing seams, exactly as on the regular path. + (1-D) and rectilinear (2-D) share one path. Ragged fields ride the SAME + per-object slab pass (issue #209): their vlen-bytes array shares the dense + ShardingCodec geometry, so a shard's digests collapse to one object instead + of the old ~K×7 per-inner-chunk CSR objects. Only the ``resolution: chunk`` + companions stay per inner chunk (their arrays are one block per chunk on the + coarse chunk grid, never sharded). Parameters ---------- @@ -327,28 +312,29 @@ def _local_region(block_index): # Group the dispatch shard's inner chunks by sharding object so each object is # accumulated and written independently; preserve encounter order for stable, - # deterministic writes. Companions + ragged stay per inner chunk (unsharded). + # deterministic writes. Companions stay per inner chunk (unsharded); ragged + # fields ride the SAME per-object slab pass as the dense arrays (issue #209 — + # their vlen-bytes array shares the dense ShardingCodec geometry), so the + # per-inner-chunk CSR fan-out (issues #142/#186) is gone. objects: dict = {} - ragged_writes: list = [] for block_index, carrier, ragged in chunk_results: - if not _carrier_empty(carrier): - objects.setdefault(_object_block(block_index), []).append((block_index, carrier)) + if not _carrier_empty(carrier) or ragged: + objects.setdefault(_object_block(block_index), []).append( + (block_index, carrier, ragged) + ) # Companions for this inner chunk are not sharded: write them straight - # through (one block per chunk). Ragged (CSR) subgroups are collected and - # fanned out below (issue #142) -- they target disjoint prefixes, so the - # per-chunk serial loop needlessly dominated a t-digest shard's write time. + # through (one block per chunk) — a chunk-resolution ragged companion + # included (its vlen array is one block per chunk on the chunk grid). if chunk_res_fields: _write_companion_columns(carrier, store, grid, block_index, chunk_res_fields) - if ragged: - ragged_writes.append((ragged, _block_index_key(block_index, grid))) - - _write_ragged_fanout(ragged_writes, store, grid=grid) + _write_ragged_companions(ragged, store, grid, block_index, chunk_res_fields) # One accumulate→write→free pass per sharding object: each holds at most one # object's slab resident at a time (the phase-8 memory bound). for obj_block, members in objects.items(): slabs: dict = {} - for block_index, carrier in members: + ragged_slabs: dict = {} + for block_index, carrier, ragged in members: region = _local_region(block_index) for name, values in _iter_carrier_columns(carrier): if name in chunk_res_fields: @@ -365,6 +351,13 @@ def _local_region(block_index): fill = array.metadata.fill_value slabs[name] = np.full((*slab_shape, *trailing), fill, dtype=values.dtype) slabs[name][region] = values.reshape((*inner_shape, *trailing)) + # Cell-resolution ragged fields: encode this chunk's payloads into + # the object slab at the same region the dense columns use. Empty + # cells keep the b"" fill, so an all-empty inner chunk is omitted + # from the shard index (sub-shard sparsity, issue #209). + _accumulate_ragged_slabs( + ragged, ragged_slabs, region, grid, slab_shape, chunk_res_fields + ) # One block selection per dense array == one shard object write. for name, slab in slabs.items(): @@ -390,79 +383,15 @@ def _local_region(block_index): f"{trailing} (the payload dim must be one whole chunk)" ) array.set_block_selection(block_idx, slab) - return store - -def _write_ragged_fanout(ragged_writes: list, store: Store, *, grid) -> None: - """Write the sharded shard's K ragged (CSR) subgroups concurrently (issue #142). - - On the sharded path :func:`write_shard_to_zarr` already holds all K inner-chunk - carriers resident (it bundles the dense side into one shard object), and each - inner chunk emits an independent CSR subgroup at a disjoint prefix - (``{group_path}/{field}/{block_key}/...``). Those writes are therefore - embarrassingly parallel and latency-bound (~K×3 tiny objects, tens of MB), so a - serial loop dominated a t-digest shard's write time. Fan them out over a bounded - ``ThreadPoolExecutor``; the on-disk layout is unchanged (still one CSR group per - inner chunk -- only the write scheduling differs), and concurrent writes to - disjoint keys of one store are already zagg's model (the dispatcher fans cells - over threads onto a shared store). - - The bound is ``min(_RAGGED_WRITE_CONCURRENCY, len(writes), fd_safe_max_workers())`` - so this per-worker inner fan-out respects the same open-file ceiling the Lambda - dispatcher guards (:func:`zagg.concurrency.fd_safe_max_workers`) -- each in-flight - write holds a store socket. A failure in any subgroup is surfaced (not swallowed): - the first exception is re-raised after the pool drains, tagged with how many of - the K writes failed. - - (The non-sharded/streaming path deliberately stays serial: it writes-then-frees - each chunk to bound peak memory (issue #91), so it never holds the K carriers at - once to parallelize.) - """ - if not ragged_writes: - return - workers = min(_RAGGED_WRITE_CONCURRENCY, len(ragged_writes), fd_safe_max_workers()) - if workers <= 1: - for ragged, ragged_key in ragged_writes: - write_ragged_to_zarr(ragged, store, grid=grid, shard_key=ragged_key) - return - - errors: list = [] - with ThreadPoolExecutor(max_workers=workers) as pool: - futures = { - pool.submit(write_ragged_to_zarr, ragged, store, grid=grid, shard_key=ragged_key): ( - ragged_key - ) - for ragged, ragged_key in ragged_writes - } - for fut in as_completed(futures): - exc = fut.exception() - if exc is not None: - errors.append((futures[fut], exc)) - if errors: - first_key, first_exc = errors[0] - # Observability (issue #186 ask 1): the chained cause never reaches - # the operator — the worker envelope and CloudWatch both record only - # this summary string — so log the first failure's full traceback - # here, and carry its type + message in the summary text too. A - # bounded one-line-per-failure roll call (after the header, so the - # indented lines read as its continuation) makes the race signature - # ACROSS subgroups visible (which keys, same or different errors). - logger.error( - f"ragged (CSR) subgroup write failed (shard_key {first_key}): " - f"{type(first_exc).__name__}: {first_exc}", - exc_info=first_exc, - ) - if len(errors) > 1: - for key, exc in errors[:_RAGGED_FAILURE_LOG_CAP]: - logger.error(f" ragged (CSR) subgroup {key}: {type(exc).__name__}: {exc}") - if len(errors) > _RAGGED_FAILURE_LOG_CAP: - n_more = len(errors) - _RAGGED_FAILURE_LOG_CAP - logger.error(f" ... and {n_more} more subgroup write failure(s)") - raise RuntimeError( - f"ragged (CSR) write failed for {len(errors)} of {len(ragged_writes)} " - f"subgroup(s) on the sharded path (first failing shard_key {first_key}: " - f"{type(first_exc).__name__}: {first_exc})" - ) from first_exc + # ONE object write per ragged field too (issue #209): the vlen-bytes + # array shares the dense arrays' ShardingCodec geometry, so the whole + # slab lands at the same object block — the ShardingCodec emits a + # single object with its internal index in place of the old ~K×7 + # per-inner-chunk CSR objects. + for name, slab in ragged_slabs.items(): + _set_ragged_block(store, grid, name, tuple(obj_block), slab) + return store def _write_companion_columns(carrier, store, grid, block_index, chunk_res_fields): @@ -496,139 +425,245 @@ def write_ragged_to_zarr( store: Store, *, grid, - shard_key: int | str, + chunk_idx: tuple, ) -> Store: - """Write a shard's ``kind: ragged`` (CSR) fields to the Zarr store (issue #48). - - Mirrors the :func:`write_dataframe_to_zarr` seam for the dense path: the - worker collects the per-cell variable-length payloads (a ``kind: ragged`` - field has no fixed per-cell width, so it cannot ride the dense block writer), - and this function persists them via :func:`zagg.csr.write_csr` — one CSR group - (``values`` / ``offsets`` / ``cell_ids``) per field per shard. - - Store layout (the contract the ``readers/tdigest_tensor.py`` reader consumes):: - - {group_path}/{field}/{shard_key}/values - {group_path}/{field}/{shard_key}/offsets - {group_path}/{field}/{shard_key}/cell_ids - - At **cell resolution** (default) ``cell_ids[k]`` is each populated cell's - position in the chunk's ``children`` block (the index collected by - ``process_shard``); the per-shard subgroup name is the ``shard_key`` label - (the coverage cell's **decimal morton string** for HEALPix — issue #199, D1 - in ``docs/design/sparse_coverage.md``), recovered by the reader directly - from the store. + """Write one chunk's ``kind: ragged`` fields to their vlen-bytes arrays (#209). + + The sharded vlen-bytes layout (issue #209) replaces the per-inner-chunk CSR + subgroups: each ragged field is ONE ``variable_length_bytes`` array on the + cell grid (template: ``grids.base.ragged_array_spec``), and each populated + cell's value is the raw little-endian bytes of its ``(n, *inner_shape)`` + payload — the interpretation recorded in the array's attrs + (``grids.base.RAGGED_ELEMENT_ATTR``). Empty cells keep the ``b""`` fill, so + a chunk with no ragged data leaves its inner chunk absent on disk. + + This is the per-chunk seam, mirroring :func:`write_dataframe_to_zarr`: it + writes the chunk's cells at storage block ``chunk_idx``. It is used on the + UNSHARDED per-chunk write paths (the runner / Lambda ``_write_chunk`` + streaming callback), where the ragged array is regular-chunked — one object + per inner chunk, so per-chunk writes stay independent (no read-modify-write + of a shared shard object). The sharded flat path bundles all K chunks in + :func:`write_shard_to_zarr`; the hive leaf in + :func:`write_ragged_leaf_to_zarr` — one object per shard on both. At **chunk resolution** (``resolution: chunk``, issue #82) a ragged field - stores ONE variable-length payload per chunk, not per cell. The populated - cells are collapsed to that single chunk payload under the same chunk-uniform - contract as scalar/vector chunk companions (every populated cell must carry an - identical payload, else raise); it is written as a single-entry CSR with - ``cell_ids == [0]`` (the lone chunk), so the on-disk layout is the same three - arrays — a consumer reads the chunk payload as the only populated "cell". + stores ONE payload per chunk: the populated cells collapse under the same + chunk-uniform contract as scalar/vector companions, and the single payload's + bytes land at the chunk's block of the chunk-grid companion array. Parameters ---------- ragged : dict - ``{field_name: (values_list, cell_ids)}`` as filled by ``process_shard``'s - ``ragged_out`` sink. Empty (or all-empty payloads) writes empty CSR arrays - (``write_csr`` skips empties), so a shard with no ragged data is a clean - no-op rather than a special case. A located field (issue #87) arrives as - ``(values_list, cell_ids, locations_list)`` and additionally writes a - ``{field}/{shard_key}/locations`` uint64 array sharing the offsets. + ``{field_name: (values_list, cell_ids)}`` as filled by ``process_shard``; + ``cell_ids`` are positions in THIS chunk's ``children`` block. A located + field (issue #87) arrives as ``(values_list, cell_ids, locations_list)`` + and additionally writes the sibling ``{field}_locations`` uint64 vlen + array, row-aligned with the payload. store : Store - Zarr-compatible store. + Zarr store with the template (including the ragged arrays) present. grid : OutputGrid - Provides ``group_path`` for routing the write (and ``config`` for the + Provides ``group_path``/``chunk_shape`` (and ``config`` for the per-field dtype + resolution). - shard_key : int or str - CSR subgroup name, used verbatim: the grid's shard label at cell - resolution (``grid.shard_label`` — decimal morton string for HEALPix, - issue #199), or the flattened block-index int at K>1. + chunk_idx : tuple of int + The chunk's storage block index (same value handed to + :func:`write_dataframe_to_zarr` for the dense columns). Returns ------- Store - The same store, with the ragged CSR arrays written. + The same store, with the chunk's ragged payloads written. """ if not ragged: return store - agg_fields = get_agg_fields(grid.config) if getattr(grid, "config", None) else {} chunk_res_fields = _chunk_resolution_fields(getattr(grid, "config", None)) - shard_key = str(shard_key) - for name, entry in ragged.items(): - # Located fields (issue #87) deliver (values_list, cell_ids, locations_list); - # unlocated fields keep the 2-tuple. - if len(entry) == 3: - values_list, cell_ids, locations_list = entry - else: - values_list, cell_ids = entry - locations_list = None - sig = get_output_signature(agg_fields[name]) if name in agg_fields else {} - dtype = sig.get("dtype") or "float32" + chunk_idx = tuple(int(i) for i in chunk_idx) + inner_shape = tuple(int(s) for s in grid.chunk_shape) + _write_ragged_companions(ragged, store, grid, chunk_idx, chunk_res_fields) + slabs: dict = {} + _accumulate_ragged_slabs( + ragged, slabs, tuple(slice(0, s) for s in inner_shape), grid, inner_shape, chunk_res_fields + ) + for name, slab in slabs.items(): + _set_ragged_block(store, grid, name, chunk_idx, slab) + return store + + +def write_ragged_leaf_to_zarr(ragged_chunks: list, store: Store, *, grid) -> Store: + """Write a hive leaf's ragged fields in ONE array write each (issue #209). + + The hive counterpart of the sharded path's slab pass: ``ragged_chunks`` is + ``[(local_block_index, ragged), ...]`` — one entry per streamed chunk, at + leaf-LOCAL blocks (``hive.leaf_block_index``). The leaf template shards a + ragged field's vlen array across the whole shard (``shard_spec``), so the K + chunks accumulate into one shard-wide object slab written in a single call + — the ShardingCodec emits ONE object per leaf in place of the per-chunk CSR + subgroups (~7 objects per populated inner chunk). A per-chunk write here + would read-modify-write that shared object K times, which is why the hive + write path collects the ragged payloads instead of streaming them — the + memory bound of that accumulation is quantified and accepted at the + collection site (``hive.process_and_write_hive``, ~200 MB peak at the o8 + t-digest scale). + + ``resolution: chunk`` ragged companions are written per chunk block (their + array is one block per chunk, unsharded — same as the scalar/vector + companions). + """ + if not any(ragged for _block, ragged in ragged_chunks): + return store + chunk_res_fields = _chunk_resolution_fields(getattr(grid, "config", None)) + slab_shape = tuple(int(s) for s in grid.shard_slab_shape()) + inner_shape = tuple(int(s) for s in grid.chunk_shape) + slabs: dict = {} + for block, ragged in ragged_chunks: + block = tuple(int(b) for b in block) + _write_ragged_companions(ragged, store, grid, block, chunk_res_fields) + region = tuple(slice(b * c, (b + 1) * c) for b, c in zip(block, inner_shape)) + _accumulate_ragged_slabs(ragged, slabs, region, grid, slab_shape, chunk_res_fields) + for name, slab in slabs.items(): + _set_ragged_block(store, grid, name, (0,) * len(slab_shape), slab) + return store + + +def _ragged_entry(entry) -> tuple: + """Normalize a ragged sink entry to ``(values_list, cell_ids, locations_list)``. + + Located fields (issue #87) deliver the 3-tuple; unlocated fields keep the + 2-tuple contract (``locations_list`` is ``None``). + """ + if len(entry) == 3: + return entry + values_list, cell_ids = entry + return values_list, cell_ids, None + + +def _ragged_sig(name: str, grid) -> dict: + """Output signature of a ragged field, tolerant of config-less stub grids.""" + agg_fields = get_agg_fields(grid.config) if getattr(grid, "config", None) else {} + if name in agg_fields: + return get_output_signature(agg_fields[name]) + return {"kind": "ragged", "inner_shape": (), "dtype": None} + + +def _ragged_payload_bytes(name: str, value, sig: dict) -> bytes: + """Raw little-endian bytes of one cell's ragged payload (the vlen wire value). + + The payload must tile the field's declared ``inner_shape`` (the guard + ``write_csr`` enforced structurally); the byte order is pinned little-endian + so the stored value matches the ``RAGGED_ELEMENT_ATTR`` interpretation on + any producer. + """ + dtype = np.dtype(sig.get("dtype") or "float32").newbyteorder("<") + arr = np.asarray(value, dtype=dtype) + inner = sig.get("inner_shape") or () + width = int(np.prod(inner)) if inner else 1 + if arr.size % width: + raise ValueError( + f"ragged field {name!r}: cell payload of {arr.size} elements does not " + f"tile inner_shape {tuple(inner)}" + ) + return np.ascontiguousarray(arr).tobytes() + + +def _accumulate_ragged_slabs(ragged, slabs, region, grid, slab_shape, chunk_res_fields) -> None: + """Encode one chunk's cell-resolution ragged payloads into the object slabs. + + ``slabs`` maps array name → object slab of ``slab_shape`` (created lazily, + ``b""``-filled — the vlen fill, so untouched cells stay absent). The chunk's + cells land in ``slabs[...][region]`` at their ``cell_ids`` position within + the chunk (row-major over ``grid.chunk_shape``, grid-agnostic via + ``np.unravel_index``). A located field (issue #87) fills the sibling + ``{name}_locations`` slab row-aligned with the payload, with the same + per-cell length contract ``write_csr`` enforced. + """ + inner_shape = tuple(int(s) for s in grid.chunk_shape) + for name, entry in (ragged or {}).items(): if name in chunk_res_fields: - # resolution: chunk — collapse the populated cells to the single chunk - # payload (chunk-uniform, like scalar/vector companions) and store it as - # a one-entry CSR (the lone chunk at cell_ids == [0]). - # ``location`` + ``resolution: chunk`` is rejected at config validation, - # but a config built without validate_config (direct PipelineConfig / - # Lambda dict payload) could still deliver a located triple here — fail - # loudly rather than silently dropping the location channel (issue #87). - if locations_list is not None: + continue # chunk companion — written per chunk block + values_list, cell_ids, locations_list = _ragged_entry(entry) + if len(values_list) != len(cell_ids): + raise ValueError( + f"values_list (len {len(values_list)}) and cell_ids (len {len(cell_ids)}) " + "must have the same length" + ) + sig = _ragged_sig(name, grid) + if name not in slabs: + slabs[name] = np.full(slab_shape, b"", dtype=object) + out = slabs[name][region] + loc_out = None + if locations_list is not None: + if len(locations_list) != len(values_list): raise ValueError( - f"ragged field {name!r} is resolution: chunk but carries a location " - f"channel; located ragged fields are cell-resolution only" + f"locations_list (len {len(locations_list)}) and values_list " + f"(len {len(values_list)}) must have the same length" ) - chunk_payload = _chunk_uniform_ragged(name, values_list) - if chunk_payload is None: - continue # whole chunk is fill — nothing to record - write_csr( - store, - f"{grid.group_path}/{name}/{shard_key}", - [chunk_payload], - [0], - dtype=dtype, - ) + loc_name = ragged_locations_name(name) + if loc_name not in slabs: + slabs[loc_name] = np.full(slab_shape, b"", dtype=object) + loc_out = slabs[loc_name][region] + loc_sig = {"kind": "ragged", "inner_shape": (), "dtype": "uint64"} + width = int(np.prod(sig.get("inner_shape") or ())) # prod(()) == 1 + locs = locations_list if locations_list is not None else [None] * len(values_list) + for cid, value, loc in zip(cell_ids, values_list, locs): + payload = _ragged_payload_bytes(name, value, sig) + pos = np.unravel_index(int(cid), inner_shape) + if loc is not None and loc_out is not None: + n_rows = np.asarray(value).size // width + loc_arr = np.asarray(loc) + if loc_arr.shape != (n_rows,): + raise ValueError( + f"locations for cell {int(cid)} of {name!r} have shape " + f"{loc_arr.shape}, expected ({n_rows},) to stay row-aligned " + f"with the payload" + ) + if n_rows: + loc_out[pos] = _ragged_payload_bytes(loc_name, loc_arr, loc_sig) + if payload: + out[pos] = payload + + +def _write_ragged_companions(ragged, store, grid, block_index, chunk_res_fields) -> None: + """Write a chunk's ``resolution: chunk`` RAGGED companions (issue #82/#209). + + The vlen analogue of :func:`_write_companion_columns`: the populated cells + collapse to the single chunk payload (chunk-uniform contract, + :func:`_chunk_uniform_ragged`), whose bytes land at the chunk's block of + the chunk-grid vlen array. ``location`` + ``resolution: chunk`` is rejected + at config validation, but a config built without ``validate_config`` could + still deliver a located triple here — fail loudly rather than silently + dropping the location channel (issue #87). + """ + if not chunk_res_fields: + return + block_idx = tuple(int(i) for i in block_index) + for name, entry in (ragged or {}).items(): + if name not in chunk_res_fields: continue - write_csr( - store, - f"{grid.group_path}/{name}/{shard_key}", - values_list, - cell_ids, - dtype=dtype, - locations_list=locations_list, + values_list, _cell_ids, locations_list = _ragged_entry(entry) + if locations_list is not None: + raise ValueError( + f"ragged field {name!r} is resolution: chunk but carries a location " + f"channel; located ragged fields are cell-resolution only" + ) + chunk_payload = _chunk_uniform_ragged(name, values_list) + if chunk_payload is None: + continue # whole chunk is fill — nothing to record + block = np.full((1,) * len(block_idx), b"", dtype=object) + block[(0,) * len(block_idx)] = _ragged_payload_bytes( + name, chunk_payload, _ragged_sig(name, grid) ) - return store - + _set_ragged_block(store, grid, name, block_idx, block) -def _block_index_key(block_index, grid) -> int: - """Flatten a chunk's block-index tuple to the CSR subgroup key (issue #48, K>1). - 1-D grids (the HEALPix companion grid, the typical case) yield a single-element - block index used directly; a multi-axis (rectilinear) block index is packed - row-major against the grid's ``chunk_grid_shape`` so each chunk maps to a - distinct CSR subgroup name. Deriving the per-axis strides from the chunk grid - (rather than a fixed shift) keeps the pack injective for any grid size. - - Shared by both the local runner and the Lambda handler (issue #82 phase 7) so - the K>1 ragged-write key is computed identically off-Lambda and on-Lambda. - """ - block = tuple(int(b) for b in block_index) - if len(block) == 1: - return block[0] - # Row-major flatten with each axis's true extent as the stride. - shape = tuple(int(s) for s in getattr(grid, "chunk_grid_shape", ())) - if len(shape) != len(block): - # Fall back to a generous fixed stride if the grid does not expose a - # matching chunk_grid_shape (keeps a unique-enough key without crashing). - key = 0 - for b in block: - key = key * (1 << 32) + b - return key - key = 0 - for b, extent in zip(block, shape): - key = key * extent + b - return key +def _set_ragged_block(store, grid, name, block_idx, block) -> None: + """One block-selection write into a ragged vlen array (issue #209).""" + with config.set({"async.concurrency": 128}): + array = open_array( + store, + path=f"{grid.group_path}/{name}", + zarr_format=3, + consolidated=False, + ) + array.set_block_selection(tuple(block_idx), block) def _chunk_uniform_ragged(name: str, values_list: list): diff --git a/src/zagg/readers/__init__.py b/src/zagg/readers/__init__.py index c537e01..f75f7ff 100644 --- a/src/zagg/readers/__init__.py +++ b/src/zagg/readers/__init__.py @@ -10,9 +10,17 @@ from zagg.readers.tdigest_tensor import ( chunk_z_range, rasterize_cell, + read_cell, read_locations, read_raw_values, read_tensors, ) -__all__ = ["chunk_z_range", "rasterize_cell", "read_locations", "read_raw_values", "read_tensors"] +__all__ = [ + "chunk_z_range", + "rasterize_cell", + "read_cell", + "read_locations", + "read_raw_values", + "read_tensors", +] diff --git a/src/zagg/readers/tdigest_tensor.py b/src/zagg/readers/tdigest_tensor.py index 4f8dee3..12374d2 100644 --- a/src/zagg/readers/tdigest_tensor.py +++ b/src/zagg/readers/tdigest_tensor.py @@ -1,37 +1,58 @@ -"""Client-side reader: per-cell t-digests → fixed ``(64, 64, n_bins)`` tensors. +"""Client-side reader: per-cell t-digests → fixed ``(side, side, n_bins)`` tensors. These are *read helpers* (issue #79) that live in-repo but sit just outside zagg's core write path: a client consuming a gridded Zarr product wants a dense, fixed-size tensor per coverage chunk, reconstructed from the per-cell t-digest ragged field that zagg writes (issue #48). -Store layout consumed ---------------------- -The t-digest field is stored per shard as a CSR ragged group (see -:mod:`zagg.csr`): under the field prefix, one subgroup per shard, named by the -shard's **parent morton id** (the coverage cell), each holding the three CSR -arrays ``values`` / ``offsets`` / ``cell_ids``:: - - {field}/{parent_morton}/values - {field}/{parent_morton}/offsets - {field}/{parent_morton}/cell_ids - -``cell_ids[k]`` is a cell's position in the chunk's row-major ``(64, 64)`` -children block, and ``values[offsets[k]:offsets[k+1]]`` is that cell's -``(k_centroids, 2)`` ``(mean, weight)`` digest. The subgroup name is the -chunk's morton index — its **decimal morton string** since issue #199 (D1 in -``docs/design/sparse_coverage.md``) — recovered directly from the store (issue -#79 design decision (3)); when a sibling ``{field}/morton`` ``uint64`` -coordinate array is present it is used in preference, mapping chunk order → -morton id. +Store layout consumed (issue #209 — the sharded vlen-bytes layout) +------------------------------------------------------------------ +A ragged field is ONE ``variable_length_bytes`` array on the cell grid:: + + {group}/{field} <- vlen array; cell i holds the raw little-endian + bytes of its (n, *inner_shape) payload + {group}/{field}_locations <- located fields only (issue #87): the uint64 + per-row location words, row-aligned + {group}/morton <- per-cell uint64 morton coordinate (zagg's + standard HEALPix coordinate array) + +The element interpretation is self-describing via the array attrs +(``grids.base.RAGGED_ELEMENT_ATTR``):: + + attrs["ragged"] = {"element": {"dtype": "float32", "shape": [-1, 2]}, + "locations": ""} # located fields only + +so the readers decode what the writer declared rather than hardcoding a dtype, +and bind the location channel by metadata, not naming convention. A store +without these attrs is not a zagg ragged vlen array (pre-issue-209 CSR stores +are a hard break) and raises a pointed error. + +The read plan honors the layout the writer chose: a whole-store sweep LISTs the +array's stored chunk objects (one object per shard under the ShardingCodec, or +one per inner chunk on the unsharded flat layout — both self-describing in the +array metadata, so one code path reads either), then decodes per read chunk. +Each read chunk is a square ``(side, side)`` block of cells (row-major +``cell_id`` within the chunk, ``side = isqrt(cells_per_chunk)`` — 64 for the +production ``chunk_inner`` configs); its coverage-cell morton id is derived +from the sibling ``morton`` coordinate. Random access to one cell +(:func:`read_cell`) indexes the vlen array directly — 2 ranged GETs on a +sharded store (index suffix + one inner chunk), never the whole shard. + +**Hive products are read one leaf at a time**: a leaf zarr (issue #199) is +exactly this layout scoped to one shard — the same ``{group}`` path, the +``morton`` sibling, the versioned ragged attrs, and the whole-leaf +ShardingCodec (one stored span) — so open the leaf store +(``hive.shard_leaf_path``) and pass the same ``field`` path. The readers are +store-scoped and never traverse the hive digit tree (leaf discovery is the +walker's/coverage MOC's job, issue #200). The reader (generator) ---------------------- -:func:`read_tensors` opens the store and yields ``(tensor, morton_index)`` one -chunk at a time. Per chunk it derives a tail-trimmed z-range from the cells' -``bottom``/``top`` quantiles, anchors a fixed ``n_bins * resolution`` window at -the floor of that range, and rasterizes each cell's digest into per-bin counts -via :func:`zagg.stats.tdigest.cdf_from_tdigest`. +:func:`read_tensors` yields ``(tensor, morton_index)`` one chunk at a time. +Per chunk it derives a tail-trimmed z-range from the cells' ``bottom``/``top`` +quantiles, anchors a fixed ``n_bins * resolution`` window at the floor of that +range, and rasterizes each cell's digest into per-bin counts via +:func:`zagg.stats.tdigest.cdf_from_tdigest`. """ from __future__ import annotations @@ -44,8 +65,7 @@ import zarr from zarr.abc.store import Store -from zagg.csr import iter_csr_cells, read_csr -from zagg.grids.morton import morton_word +from zagg.grids.base import RAGGED_ELEMENT_ATTR, RAGGED_SPEC from zagg.stats.tdigest import cdf_from_tdigest, quantile_from_tdigest __all__ = [ @@ -54,12 +74,9 @@ "read_tensors", "read_raw_values", "read_locations", + "read_cell", ] -# Coverage chunk is a 64×64 block of child cells (issue #79). -_CHUNK_SIDE = 64 -_CHUNK_CELLS = _CHUNK_SIDE * _CHUNK_SIDE - FitMode = Literal["raise", "degrade_resolution", "collapse_bins"] TensorDtype = Literal["uint16", "uint32", "float32"] @@ -212,86 +229,163 @@ def chunk_z_range( raise ValueError(f"unknown fit mode {fit!r}") -def _is_decimal_shard_name(name: str) -> bool: - """Whether a CSR subgroup name fits the decimal morton *shard-id* grammar. - - Grammar: optional sign, base digit ``1..6``, then one ``1..4`` digit per - order — with at least ONE order digit. zagg shard orders are >= 1, so a - real shard id always has >= 2 digits; that floor is what disambiguates the - lone digits ``"1".."6"``, which the raw grammar would accept as order-0 ids - but which in practice are small block-index keys (review finding, PR #205). - """ - body = name[1:] if name.startswith("-") else name - return len(body) >= 2 and body[:1] in "123456" and all(d in "1234" for d in body[1:]) +# --------------------------------------------------------------------------- # +# vlen-bytes store access (issue #209) +# --------------------------------------------------------------------------- # -def _subgroup_words(shard_keys: list[str]) -> dict[str, int]: - """Parse CSR subgroup names to their numeric keys (issue #199). +def _open_ragged(store: Store, field: str, zarr_format) -> tuple: + """Open a ragged vlen array; return ``(arr, elem_dtype, elem_shape, meta)``. - Shard-keyed subgroups are named by the decimal morton string (D1), so the - numeric key is the parsed packed word; block-index-keyed subgroups (the K>1 - layout) keep plain non-negative ints. The flavor is decided per store: - - - every name fits the shard-id grammar (:func:`_is_decimal_shard_name`) - -> decimal morton store; keys are the parsed packed words; - - otherwise -> block-keyed store; keys are ``int(name)`` — and a leading - ``-`` in this branch raises, because a signed name is provably NOT a - block index (block keys are non-negative flattened indices): a mixed - store (e.g. pre-#199 debris alongside decimal-named subgroups) must fail - loudly rather than silently mis-key every shard. + ``meta`` is the array's :data:`~zagg.grids.base.RAGGED_ELEMENT_ATTR` attrs + payload — the element interpretation the WRITER declared, which is what the + readers decode by (never a hardcoded dtype). Raises a pointed error when + the attrs are missing/malformed: silently guessing an element layout would + misinterpret every payload (pre-issue-209 CSR stores are a hard break). """ - if shard_keys and all(_is_decimal_shard_name(k) for k in shard_keys): - return {k: morton_word(k) for k in shard_keys} - signed = [k for k in shard_keys if k.startswith("-")] - if signed: + arr = zarr.open_array(store, path=field, mode="r", zarr_format=zarr_format) + raw_meta = arr.attrs.get(RAGGED_ELEMENT_ATTR) + meta: dict = dict(raw_meta) if isinstance(raw_meta, dict) else {} + element = meta.get("element") + if not isinstance(element, dict) or "dtype" not in element or "shape" not in element: + raise ValueError( + f"{field!r} carries no ragged element declaration " + f'(attrs["{RAGGED_ELEMENT_ATTR}"]["element"]); it is not a zagg ragged ' + f"vlen array (issue #209). Pre-issue-209 CSR stores are not readable " + f"— rewrite the store." + ) + # Version gate (the coverage-envelope discipline): an unknown/future spec + # must fail loudly, never half-parse. This attrs seam is the INTERIM + # contract the issue #210 typed vlen-array dtype migration supersedes. + if meta.get("spec") != RAGGED_SPEC: + raise ValueError( + f"{field!r} declares ragged spec {meta.get('spec')!r}; this reader " + f"understands {RAGGED_SPEC!r} only — a newer writer's layout must be " + f"adopted deliberately, not half-parsed" + ) + if arr.ndim != 1: raise ValueError( - f"CSR subgroup names mix decimal morton ids with non-id names " - f"(signed name(s) {signed[:4]} cannot be block-index keys); the " - f"store layout is ambiguous — likely pre-issue-199 debris next to " - f"decimal-named subgroups. Rewrite the store." + f"{field!r} has {arr.ndim} dimensions; the t-digest tensor readers " + f"consume HEALPix products (1-D cells axis)" ) - return {k: int(k) for k in shard_keys} + try: + dtype = np.dtype(str(element["dtype"])).newbyteorder("<") + except TypeError as e: + raise ValueError( + f"{field!r} declares an unreadable element dtype " + f'{element["dtype"]!r} in attrs["{RAGGED_ELEMENT_ATTR}"]["element"]' + ) from e + shape = tuple(int(s) for s in element["shape"]) + return arr, dtype, shape, meta -def _resolve_chunk_morton( - store: Store, field: str, shard_keys: list[str], zarr_format: Literal[2, 3] -) -> dict[str, int]: - """Map each shard subgroup name → its morton id. +def _decode_cell(raw, dtype: np.dtype, shape: tuple) -> np.ndarray: + """One cell's payload from its raw vlen bytes, per the declared element.""" + return np.frombuffer(bytes(raw), dtype=dtype).reshape(shape) - Prefers a sibling ``{field}/morton`` ``uint64`` coordinate array (chunk - order → morton id); falls back to parsing the subgroup name as the parent - morton id (the shard key is the parent morton — see :func:`process_shard`), - a decimal morton string since issue #199. - The coordinate array is aligned against the subgroup names sorted - **numerically** by parsed key (the canonical ascending-morton chunk order), - not lexicographically — string sorting would mis-align names of differing - digit counts (e.g. ``"1000"`` before ``"99"``). +def _stored_chunk_spans(arr) -> list[tuple[int, int]]: + """Cell spans ``(start, stop)`` of the array's STORED objects, ascending. + + The whole-store read plan: one LIST of the array's ``c/`` data + keys instead of probing every chunk of a mostly-fill array. Under the + ShardingCodec an object spans a whole shard (``arr.shards``); on the + unsharded flat layout it is one read chunk — both derive from the array's + own metadata, so the sharded and per-inner-chunk layouts share this path. + ``zarr.core.sync.sync`` runs the store's async listing on zarr's own event + loop (zarr 3 exposes no public sync listing). + """ + from zarr.core.sync import sync + + async def _collect(gen): + return [key async for key in gen] + + span = int((arr.shards or arr.chunks)[0]) + prefix = f"{arr.path}/c/" if arr.path else "c/" + keys = sync(_collect(arr.store_path.store.list_prefix(prefix))) + ordinals = sorted(int(k.rsplit("/", 1)[-1]) for k in keys) + return [(o * span, min((o + 1) * span, int(arr.shape[0]))) for o in ordinals] + + +def _iter_populated_chunks(arr) -> Iterator[tuple[int, list]]: + """Yield ``(chunk_start, [(cell_pos, raw_bytes), ...])`` per populated chunk. + + Restricted to the stored objects (:func:`_stored_chunk_spans`), each read + in ONE slice: slicing per inner chunk would re-fetch a sharded object's + index suffix on every ``__getitem__`` (~K redundant GETs per shard at the + production K=256 — review, PR #211), while the full-span slice fetches + the stored object ONCE (zarr reads a whole outer chunk in a single GET + and splits it locally). The held cost is one span's decoded payload — + ~141 MB at the o8 t-digest scale, the same bound the hive write side + documents and accepts (``process_and_write_hive``), and this is the + client-side bulk reader. Chunks whose cells are all absent (the ``b""`` + fill) are skipped; ``cell_pos`` is the cell's row-major position within + the chunk — the same index the writer placed it at. """ - words = _subgroup_words(shard_keys) + cells_per_chunk = int(arr.chunks[0]) + for span_start, span_stop in _stored_chunk_spans(arr): + span = arr[span_start:span_stop] + for offset in range(0, span_stop - span_start, cells_per_chunk): + block = span[offset : offset + cells_per_chunk] + populated = [(pos, block[pos]) for pos in range(len(block)) if len(block[pos])] + if populated: + yield span_start + offset, populated + + +def _open_morton(store: Store, field: str, zarr_format): + """The sibling per-cell ``morton`` coordinate array (chunk identity source).""" + parent, _, _name = field.rpartition("/") + path = f"{parent}/morton" if parent else "morton" try: - arr = zarr.open_array(store, path=f"{field}/morton", mode="r", zarr_format=zarr_format) - morton = np.asarray(arr[...]) - except (FileNotFoundError, KeyError, ValueError): - return words - if len(morton) != len(shard_keys): - # Coordinate present but not 1:1 with subgroups — fall back to the names. - return words - return {k: int(m) for k, m in zip(sorted(shard_keys, key=words.__getitem__), morton)} - - -def _shard_groups( - store: Store, field: str, zarr_format: Literal[2, 3] -) -> tuple[list[str], dict[str, int]]: - """Enumerate a field's per-shard CSR subgroups and their morton ids. - - The shared reader preamble: subgroup names under ``field`` (the ``morton`` - coordinate array is a sibling, not a shard) plus the name → morton map from - :func:`_resolve_chunk_morton`. + return zarr.open_array(store, path=path, mode="r", zarr_format=zarr_format) + except (FileNotFoundError, KeyError) as e: + raise ValueError( + f"{field!r} has no sibling 'morton' coordinate array at {path!r}; the " + f"vlen readers derive each chunk's coverage-cell id from the per-cell " + f"morton coordinate (issue #209)" + ) from e + + +def _chunk_word(words: np.ndarray, field: str, start: int) -> int: + """A read chunk's coverage-cell morton id from its cells' morton words. + + ``words`` are the chunk's per-cell morton coordinates (packed uint64 area + words at the cell order; ``0`` is the unwritten fill). The chunk id is any + written cell's word coarsened to the chunk order — ``cell_order - + log4(cells_per_chunk)`` — the same parent cell the CSR layout named its + subgroups by. """ - group = zarr.open_group(store, path=field, mode="r", zarr_format=zarr_format) - shard_keys = sorted(k for k in group.group_keys() if k != "morton") - return shard_keys, _resolve_chunk_morton(store, field, shard_keys, zarr_format) + from mortie import clip2order + + from zagg.grids.morton import morton_decimal + + written = words[words != 0] + if written.size == 0: + raise ValueError( + f"chunk at cell {start} of {field!r} has ragged payloads but no written " + f"'morton' coordinate — the dense coordinate write did not cover it" + ) + decimal = morton_decimal(int(written[0])) + cell_order = len(decimal) - (2 if decimal.startswith("-") else 1) + depth = (int(len(words)).bit_length() - 1) // 2 # log4(cells_per_chunk) + return int(clip2order(cell_order - depth, written[:1])[0]) + + +def _tensor_side(arr, field: str) -> int: + """Row-major square side of one read chunk (64 for the production configs).""" + cells_per_chunk = int(arr.chunks[0]) + side = math.isqrt(cells_per_chunk) + if side * side != cells_per_chunk: + raise ValueError( + f"{field!r} read chunk holds {cells_per_chunk} cells — not a square " + f"block, so it cannot rasterize to a (side, side, n_bins) tensor" + ) + return side + + +# --------------------------------------------------------------------------- # +# public readers +# --------------------------------------------------------------------------- # def read_tensors( @@ -308,18 +402,20 @@ def read_tensors( ) -> Iterator[tuple[np.ndarray, int]]: """Yield ``(tensor, morton_index)`` per coverage chunk of a t-digest field. - Opens the CSR ragged store for ``field`` and iterates one Zarr chunk (one - 64×64 parent block) at a time. For each chunk it trims the per-cell tails, - derives a fixed z-window (see :func:`chunk_z_range`), rasterizes every - populated cell's digest into ``n_bins`` counts (see :func:`rasterize_cell`), - and emits the ``(64, 64, n_bins)`` tensor with the chunk's morton id. + Sweeps the field's vlen array one read chunk (one square cell block) at a + time, visiting only the STORED objects. For each populated chunk it trims + the per-cell tails, derives a fixed z-window (see :func:`chunk_z_range`), + rasterizes every populated cell's digest into ``n_bins`` counts (see + :func:`rasterize_cell`), and emits the ``(side, side, n_bins)`` tensor with + the chunk's coverage-cell morton id (``side`` is 64 for the production + ``chunk_inner`` configs). Parameters ---------- store : Store - Zarr store holding the per-shard CSR groups under ``field``. + Zarr store holding the ragged vlen array (issue #209 layout). field : str - Field prefix (e.g. ``"h_tdigest"``). + Array path (e.g. ``"19/h_tdigest"``). n_bins : int, optional Number of z-bins (default 128). resolution : float, optional @@ -341,27 +437,29 @@ def read_tensors( Yields ------ (tensor, morton_index) : (ndarray, int) - ``tensor`` has shape ``(64, 64, n_bins_out)`` and the requested dtype; - ``morton_index`` is the chunk's coverage-cell morton id. + ``tensor`` has shape ``(side, side, n_bins_out)`` and the requested + dtype; ``morton_index`` is the chunk's coverage-cell morton id. Raises ------ ValueError - On an unknown ``dtype``/``fit``, or (with ``fit="raise"``) a chunk whose - trimmed range overflows the fixed window. + On an unknown ``dtype``/``fit``, a store missing the ragged element + attrs or the ``morton`` sibling, or (with ``fit="raise"``) a chunk + whose trimmed range overflows the fixed window. """ if dtype not in _TENSOR_DTYPES: raise ValueError(f"unknown dtype {dtype!r}; expected one of {sorted(_TENSOR_DTYPES)}") out_dtype = _TENSOR_DTYPES[dtype] is_float = np.issubdtype(out_dtype, np.floating) - shard_keys, morton_of = _shard_groups(store, field, zarr_format) + arr, elem_dtype, elem_shape, _meta = _open_ragged(store, field, zarr_format) + morton = _open_morton(store, field, zarr_format) + side = _tensor_side(arr, field) - for key in shard_keys: - cells = iter_csr_cells(read_csr(store, f"{field}/{key}", zarr_format=zarr_format)) - digests = [payload for _, payload in cells] + for start, populated in _iter_populated_chunks(arr): + cells = [(pos, _decode_cell(raw, elem_dtype, elem_shape)) for pos, raw in populated] z_lo, n_bins_c, resolution_c = chunk_z_range( - digests, + [digest for _pos, digest in cells], n_bins=n_bins, resolution=resolution, bottom=bottom, @@ -369,19 +467,15 @@ def read_tensors( fit=fit, ) - tensor = np.zeros((_CHUNK_SIDE, _CHUNK_SIDE, n_bins_c), dtype=out_dtype) + tensor = np.zeros((side, side, n_bins_c), dtype=out_dtype) for cell_id, digest in cells: - if not (0 <= cell_id < _CHUNK_CELLS): - raise ValueError( - f"cell_id {cell_id} out of range for a {_CHUNK_SIDE}×{_CHUNK_SIDE} chunk" - ) counts = rasterize_cell(digest, z_lo, resolution_c, n_bins_c) if not is_float: counts = np.rint(counts) - row, col = divmod(cell_id, _CHUNK_SIDE) + row, col = divmod(cell_id, side) tensor[row, col, :] = counts.astype(out_dtype) - yield tensor, morton_of[key] + yield tensor, _chunk_word(morton[start : start + side * side], field, start) def read_raw_values( @@ -401,9 +495,9 @@ def read_raw_values( Parameters ---------- store : Store - Zarr store holding the per-shard CSR groups under ``field``. + Zarr store holding the ragged vlen array (issue #209 layout). field : str - Field prefix. + Array path. zarr_format : int, optional Zarr format version (default 3). @@ -411,7 +505,8 @@ def read_raw_values( ------ (morton_index, cell_id, values) : (int, int, ndarray) ``values`` is the cell's recovered 1-D sample vector (sorted ascending, - as the digest stores centroids by mean). + as the digest stores centroids by mean); ``cell_id`` its row-major + position within the chunk. Raises ------ @@ -419,22 +514,21 @@ def read_raw_values( If any cell's digest carries a merged centroid (weight > 1), so the raw values cannot be recovered without loss. """ - shard_keys, morton_of = _shard_groups(store, field, zarr_format) - - for key in shard_keys: - morton = morton_of[key] - for cell_id, digest in iter_csr_cells( - read_csr(store, f"{field}/{key}", zarr_format=zarr_format) - ): - if len(digest) == 0: - continue + arr, elem_dtype, elem_shape, _meta = _open_ragged(store, field, zarr_format) + morton = _open_morton(store, field, zarr_format) + cells_per_chunk = int(arr.chunks[0]) + + for start, populated in _iter_populated_chunks(arr): + word = _chunk_word(morton[start : start + cells_per_chunk], field, start) + for cell_id, raw in populated: + digest = _decode_cell(raw, elem_dtype, elem_shape) weights = np.asarray(digest[:, 1], dtype=np.float64) if np.any(weights > 1.0): raise ValueError( - f"cell {cell_id} (chunk {morton}) has merged centroids " + f"cell {cell_id} (chunk {word}) has merged centroids " "(weight > 1); raw values are not losslessly recoverable" ) - yield int(morton), int(cell_id), np.asarray(digest[:, 0], dtype=np.float64) + yield word, int(cell_id), np.asarray(digest[:, 0], dtype=np.float64) def read_locations( @@ -446,21 +540,22 @@ def read_locations( """Yield ``(morton_index, cell_id, locations)`` per populated cell. The location-channel reader (issue #87): a located ragged field stores a - per-centroid ``uint64`` morton location vector in a ``locations`` array - sharing the CSR ``offsets``/``cell_ids`` with the field's ``values``, so - the k-th populated cell's locations are ``locations[offsets[k]:offsets[k+1]]`` - — one word per centroid, each the deepest morton cell enclosing that - centroid's member observations (an exact order-29 point word for a 1-obs - centroid). Decode or coarsen the words with mortie (``mort2healpix``, - ``clip2order``, ``common_ancestor``); geometric queries compose directly on - the packed words. + per-centroid ``uint64`` morton location vector in a sibling vlen array, + row-aligned with the digest — one word per centroid, each the deepest + morton cell enclosing that centroid's member observations (an exact + order-29 point word for a 1-obs centroid). The sibling is bound by the + payload array's attrs declaration (``attrs["ragged"]["locations"]`` — + issue #209), never by reconstructing the naming convention. Decode or + coarsen the words with mortie (``mort2healpix``, ``clip2order``, + ``common_ancestor``); geometric queries compose directly on the packed + words. Parameters ---------- store : Store - Zarr store holding the per-shard CSR groups under ``field``. + Zarr store holding the ragged vlen arrays (issue #209 layout). field : str - Field prefix. + The PAYLOAD array's path (not the sibling's). zarr_format : int, optional Zarr format version (default 3). @@ -474,33 +569,56 @@ def read_locations( Raises ------ ValueError - If the field was not written with a location channel. + If the field declares no locations channel (it was not written as a + located ragged field). """ - shard_keys, morton_of = _shard_groups(store, field, zarr_format) - - for key in shard_keys: - # Open only the three arrays this reader needs — skipping ``values`` - # halves the bytes read for a locations-only sweep (the digest payload - # is the field's largest array and is not consumed here). - prefix = f"{field}/{key}" - try: - locations = np.asarray( - zarr.open_array( - store, path=f"{prefix}/locations", mode="r", zarr_format=zarr_format - )[...] - ) - except (FileNotFoundError, KeyError) as e: - raise ValueError( - f"{prefix!r} has no locations array; it was not written as a " - f"located ragged field (declare 'location:' on the field — issue #87)" - ) from e - offsets = np.asarray( - zarr.open_array(store, path=f"{prefix}/offsets", mode="r", zarr_format=zarr_format)[...] + _arr, _dt, _sh, meta = _open_ragged(store, field, zarr_format) + sibling = meta.get("locations") + if not sibling: + raise ValueError( + f"{field!r} declares no locations channel " + f'(attrs["{RAGGED_ELEMENT_ATTR}"]["locations"]); it was not written as ' + f"a located ragged field (declare 'location:' on the field — issue #87)" ) - cell_ids = np.asarray( - zarr.open_array(store, path=f"{prefix}/cell_ids", mode="r", zarr_format=zarr_format)[ - ... - ] + parent, _, _name = field.rpartition("/") + loc_field = f"{parent}/{sibling}" if parent else str(sibling) + # Sweep the SIBLING only — skipping the payload array halves the bytes + # read for a locations-only pass (the digest payload is not consumed here). + loc_arr, loc_dtype, loc_shape, _loc_meta = _open_ragged(store, loc_field, zarr_format) + morton = _open_morton(store, field, zarr_format) + cells_per_chunk = int(loc_arr.chunks[0]) + + for start, populated in _iter_populated_chunks(loc_arr): + word = _chunk_word(morton[start : start + cells_per_chunk], loc_field, start) + for cell_id, raw in populated: + yield word, int(cell_id), _decode_cell(raw, loc_dtype, loc_shape) + + +def read_cell( + store: Store, + field: str, + cell: int, + *, + zarr_format: Literal[2, 3] = 3, +) -> np.ndarray: + """Random-access ONE cell's ragged payload, decoded per the element attrs. + + The issue #209 single-cell path: indexing the vlen array reads exactly two + ranged GETs on a sharded store — the shard-index suffix, then the one inner + chunk holding the cell — never the whole shard object. ``cell`` is the + cell's global position on the array's cells axis (no negative-index wrap). + An absent/empty cell returns the zero-length ``(0, *inner_shape)`` array; + an out-of-range index raises ``IndexError`` naming the valid range (zarr + basic selection would silently clamp the slice — review, PR #211). Works + on any zagg ragged vlen array — a located field's ``{field}_locations`` + sibling included (its elements decode as ``(n,)`` uint64 words). + """ + arr, elem_dtype, elem_shape, _meta = _open_ragged(store, field, zarr_format) + cell = int(cell) + if not 0 <= cell < int(arr.shape[0]): + raise IndexError( + f"cell {cell} out of range for {field!r} with {int(arr.shape[0])} cells " + f"(valid: 0..{int(arr.shape[0]) - 1}; negative indices do not wrap)" ) - for k, cid in enumerate(cell_ids): - yield int(morton_of[key]), int(cid), locations[offsets[k] : offsets[k + 1]] + (raw,) = arr[cell : cell + 1] + return _decode_cell(raw, elem_dtype, elem_shape) diff --git a/src/zagg/runner.py b/src/zagg/runner.py index 987f13a..a9baca8 100644 --- a/src/zagg/runner.py +++ b/src/zagg/runner.py @@ -64,7 +64,6 @@ write_ragged_to_zarr, write_shard_to_zarr, ) -from zagg.processing.write import _block_index_key from zagg.store import open_object_store, open_store logger = logging.getLogger(__name__) @@ -1099,29 +1098,22 @@ def _process_and_write( ``K = grid.chunks_per_shard`` finer Zarr chunks. ``process_shard`` reads the granules once and returns one ``(block_index, carrier, ragged)`` per chunk via ``chunk_results``; this writes each chunk's dense region (at its own - ``block_index``) plus its ragged (CSR) companion. At K==1 ``chunk_results`` has - exactly one entry whose ``block_index`` equals ``chunk_idx``, so the write is - byte-for-byte the single-chunk path. ``chunk_idx`` is retained for the K==1 - callers/signature but the per-chunk block index from ``iter_chunks`` is used. + ``block_index``) plus its ragged vlen payloads (issue #209). At K==1 + ``chunk_results`` has exactly one entry whose ``block_index`` equals + ``chunk_idx``, so the write is byte-for-byte the single-chunk path. + ``chunk_idx`` is retained for the K==1 callers/signature but the per-chunk + block index from ``iter_chunks`` is used. """ - # K==1 vs K>1 is fixed by the grid, not the materialized list (issue #91): at - # K==1 the lone chunk IS the shard so its CSR subgroup is keyed by ``shard_key``; - # at K>1 each finer chunk is keyed by its own block index. Deriving it from the - # grid lets the non-sharded path stream (no materialized count needed). - single_chunk = int(getattr(grid, "chunks_per_shard", 1)) == 1 def _write_chunk(block_index, carrier, ragged): # write_dataframe_to_zarr no-ops on an empty carrier (DataFrame or Arrow # table), so no carrier-specific emptiness check is needed here. write_dataframe_to_zarr(carrier, zarr_store, grid=grid, chunk_idx=block_index) - # Persist this chunk's ragged (CSR) fields — one CSR group per field per - # chunk (issue #48). No-ops when ``ragged`` is empty. At K==1 the CSR - # subgroup is named by the grid's shard label (the decimal morton string - # for HEALPix — issue #199); K>1 keeps the flattened block-index int. - ragged_key = ( - shard_label(grid, shard_key) if single_chunk else _block_index_key(block_index, grid) - ) - write_ragged_to_zarr(ragged, zarr_store, grid=grid, shard_key=ragged_key) + # Persist this chunk's ragged fields into their vlen-bytes arrays at the + # same block (issue #209). The array is regular-chunked on this + # unsharded path, so per-chunk writes stay independent. No-op when + # ``ragged`` is empty. + write_ragged_to_zarr(ragged, zarr_store, grid=grid, chunk_idx=block_index) # Sharded output (issue #108): the shard's K inner chunks bundle into one # ShardingCodec shard object — write the whole shard in one block selection per diff --git a/src/zagg/stats/tdigest.py b/src/zagg/stats/tdigest.py index ce01e59..753d40c 100644 --- a/src/zagg/stats/tdigest.py +++ b/src/zagg/stats/tdigest.py @@ -67,7 +67,7 @@ per-observation order-29 morton point words as ``locations=``; the reducer then returns a ``(digest, locations)`` pair whose second element is the ``(k,)`` uint64 per-centroid location — the deepest morton cell enclosing each -centroid's members (``mortie.common_ancestor``), stored as a companion CSR +centroid's members (``mortie.common_ancestor``), stored as a companion vlen array. See ``zagg/configs/atl03_tdigest_located_healpix.yaml``. """ diff --git a/tests/test_config.py b/tests/test_config.py index fa74384..f8fe36e 100644 --- a/tests/test_config.py +++ b/tests/test_config.py @@ -2016,6 +2016,49 @@ def test_location_unknown_column_rejected(self): with pytest.raises(ValueError, match="is not 'leaf_id' or a data_source variable"): validate_config(cfg) + def test_sibling_locations_field_collision_rejected(self): + # Issue #209 (review, PR #211): the located channel lands in a SIBLING + # array `{field}_locations` — a name the user never wrote. A declared + # field claiming it would silently lose in the template members dict + # and interleave into the same object slab at write time. + cfg = _ragged_cfg( + inner_shape=[2], + location="leaf_id", + function="zagg.stats.tdigest.build_tdigest", + ) + cfg.aggregation["variables"]["h_ph_tdigest_locations"] = { + "function": "min", + "source": "h_ph", + "dtype": "float32", + } + with pytest.raises(ValueError, match="sibling array named 'h_ph_tdigest_locations'"): + validate_config(cfg) + + def test_sibling_locations_coordinate_collision_rejected(self): + # Same guard for coordinates — they share the template members dict. + cfg = _ragged_cfg( + inner_shape=[2], + location="leaf_id", + function="zagg.stats.tdigest.build_tdigest", + ) + cfg.aggregation["coordinates"]["h_ph_tdigest_locations"] = { + "dtype": "uint64", + "fill_value": 0, + } + with pytest.raises(ValueError, match="sibling array named 'h_ph_tdigest_locations'"): + validate_config(cfg) + + def test_locations_suffix_free_without_location(self): + # The sibling name is reserved only for LOCATED fields: an unlocated + # ragged field coexists with a `{field}_locations`-named scalar. + cfg = _ragged_cfg(inner_shape=[2]) + cfg.aggregation["variables"]["h_ph_tdigest_locations"] = { + "function": "min", + "source": "h_ph", + "dtype": "float32", + } + validate_config(cfg) + def test_location_defaults_to_healpix_when_grid_absent(self): # No output.grid block defaults to healpix everywhere else (the grid # factory), so a located field must validate the same way (review fold). diff --git a/tests/test_csr.py b/tests/test_csr.py deleted file mode 100644 index efd6aa7..0000000 --- a/tests/test_csr.py +++ /dev/null @@ -1,205 +0,0 @@ -"""Tests for the CSR (ragged) writer/reader — issue #48, phase 3.""" - -import numpy as np -import pytest -from zarr.storage import MemoryStore - -from zagg.csr import iter_csr_cells, read_csr, write_csr - - -class TestWriteCsr: - def test_round_trip_basic(self): - """Write and read back variable-length per-cell payloads.""" - store = MemoryStore() - payloads = [ - np.array([[1.0, 2.0], [3.0, 4.0]], dtype=np.float32), # cell 0: 2 rows - np.array([[5.0, 6.0]], dtype=np.float32), # cell 2: 1 row - np.array([[7.0, 8.0], [9.0, 10.0], [11.0, 12.0]], dtype=np.float32), # cell 5: 3 rows - ] - cell_ids = [0, 2, 5] - write_csr(store, "tdigest", payloads, cell_ids) - csr = read_csr(store, "tdigest") - - assert csr["values"].shape == (6, 2) - assert len(csr["offsets"]) == 4 # n_populated + 1 - np.testing.assert_array_equal(csr["cell_ids"], [0, 2, 5]) - np.testing.assert_array_equal(csr["offsets"], [0, 2, 3, 6]) - - def test_round_trip_per_cell_slices(self): - """Per-cell reconstruction via iter_csr_cells is exact.""" - store = MemoryStore() - payloads = [ - np.array([[1.0, 10.0]], dtype=np.float32), - np.array([[2.0, 20.0], [3.0, 30.0]], dtype=np.float32), - ] - cell_ids = [7, 3] - write_csr(store, "field", payloads, cell_ids) - csr = read_csr(store, "field") - decoded = iter_csr_cells(csr) - - assert decoded[0][0] == 7 - np.testing.assert_array_almost_equal(decoded[0][1], [[1.0, 10.0]]) - assert decoded[1][0] == 3 - np.testing.assert_array_almost_equal(decoded[1][1], [[2.0, 20.0], [3.0, 30.0]]) - - def test_empty_payloads_skipped(self): - """write_csr with all-empty payloads writes empty arrays (no crash).""" - store = MemoryStore() - write_csr(store, "empty_field", [], []) - csr = read_csr(store, "empty_field") - assert csr["values"].size == 0 - np.testing.assert_array_equal(csr["offsets"], [0]) - assert csr["cell_ids"].size == 0 - - def test_mixed_empty_non_empty(self): - """Empty arrays in values_list are silently skipped.""" - store = MemoryStore() - payloads = [ - np.array([], dtype=np.float32).reshape(0, 2), - np.array([[1.0, 2.0]], dtype=np.float32), - np.array([], dtype=np.float32).reshape(0, 2), - ] - cell_ids = [0, 1, 2] - write_csr(store, "sparse", payloads, cell_ids) - csr = read_csr(store, "sparse") - # Only cell 1 has data. - np.testing.assert_array_equal(csr["cell_ids"], [1]) - np.testing.assert_array_equal(csr["offsets"], [0, 1]) - assert csr["values"].shape[0] == 1 - - def test_mismatched_lengths_raises(self): - with pytest.raises(ValueError, match="same length"): - write_csr(MemoryStore(), "f", [np.array([1.0])], [0, 1]) - - def test_inconsistent_inner_shape_raises(self): - payloads = [ - np.array([[1.0, 2.0]], dtype=np.float32), # inner (2,) - np.array([[3.0, 4.0, 5.0]], dtype=np.float32), # inner (3,) — mismatch - ] - with pytest.raises(ValueError, match="Inconsistent inner shape"): - write_csr(MemoryStore(), "f", payloads, [0, 1]) - - def test_csr_invariant_offsets_last_equals_len_values(self): - """offsets[-1] must equal len(values) (standard CSR invariant).""" - store = MemoryStore() - payloads = [np.arange(6, dtype=np.float32).reshape(3, 2)] - write_csr(store, "inv", payloads, [4]) - csr = read_csr(store, "inv") - assert csr["offsets"][-1] == len(csr["values"]) - - def test_single_cell_1d_payload(self): - """1-D per-cell arrays (scalar inner type) are written/read correctly.""" - store = MemoryStore() - payloads = [np.array([1.0, 2.0, 3.0], dtype=np.float32)] - write_csr(store, "scalar_inner", payloads, [0]) - csr = read_csr(store, "scalar_inner") - np.testing.assert_array_equal(csr["values"], [1.0, 2.0, 3.0]) - np.testing.assert_array_equal(csr["offsets"], [0, 3]) - np.testing.assert_array_equal(csr["cell_ids"], [0]) - - def test_dtype_preserved(self): - """The dtype parameter controls the values array dtype.""" - store = MemoryStore() - payloads = [np.array([[1, 2]], dtype=np.int32)] - write_csr(store, "ints", payloads, [0], dtype="int32") - csr = read_csr(store, "ints") - assert csr["values"].dtype == np.dtype("int32") - - -class TestLocationsChannel: - """Issue #87 phase 4: the uint64 ``locations`` companion CSR array shares - ``offsets``/``cell_ids`` with ``values`` and is purely additive — the - value-only three-array layout is byte-identical with and without it.""" - - _PAYLOADS = [ - np.array([[1.0, 2.0], [3.0, 4.0]], dtype=np.float32), - np.array([[5.0, 6.0]], dtype=np.float32), - ] - _CELL_IDS = [0, 5] - # Order-29-scale words (past int64 max) — uint64 carriage must be exact. - _LOCS = [ - np.array([2**63 + 11, 2**63 + 12], dtype=np.uint64), - np.array([2**61 + 5], dtype=np.uint64), - ] - - @staticmethod - def _snapshot(store): - return {k: bytes(v.to_bytes()) for k, v in store._store_dict.items()} - - def test_locations_round_trip(self): - store = MemoryStore() - write_csr(store, "f", self._PAYLOADS, self._CELL_IDS, locations_list=self._LOCS) - csr = read_csr(store, "f", locations=True) - assert csr["locations"].dtype == np.uint64 - np.testing.assert_array_equal(csr["locations"], np.concatenate(self._LOCS)) - # Shared offsets slice the locations exactly like the values. - for k in range(len(self._CELL_IDS)): - sl = slice(csr["offsets"][k], csr["offsets"][k + 1]) - np.testing.assert_array_equal(csr["locations"][sl], self._LOCS[k]) - assert len(csr["locations"][sl]) == len(csr["values"][sl]) - - def test_value_arrays_byte_identical_with_and_without_locations(self): - plain, located = MemoryStore(), MemoryStore() - write_csr(plain, "f", self._PAYLOADS, self._CELL_IDS) - write_csr(located, "f", self._PAYLOADS, self._CELL_IDS, locations_list=self._LOCS) - plain_snap = self._snapshot(plain) - located_snap = self._snapshot(located) - # The value-only store has no locations keys at all. - assert not any("locations" in k for k in plain_snap) - # The located store adds ONLY locations keys; every shared key is - # byte-identical (the regression guard for existing value-only output). - extra = set(located_snap) - set(plain_snap) - assert extra and all("locations" in k for k in extra) - for k in plain_snap: - assert located_snap[k] == plain_snap[k], f"byte drift at {k}" - - def test_empty_payloads_drop_their_locations(self): - store = MemoryStore() - payloads = [self._PAYLOADS[0], np.empty((0, 2), dtype=np.float32), self._PAYLOADS[1]] - locs = [self._LOCS[0], np.empty(0, dtype=np.uint64), self._LOCS[1]] - write_csr(store, "f", payloads, [0, 3, 5], locations_list=locs) - csr = read_csr(store, "f", locations=True) - np.testing.assert_array_equal(csr["cell_ids"], [0, 5]) - np.testing.assert_array_equal(csr["locations"], np.concatenate(self._LOCS)) - - def test_all_empty_writes_empty_locations(self): - store = MemoryStore() - write_csr( - store, - "f", - [np.empty((0, 2), dtype=np.float32)], - [0], - locations_list=[np.empty(0, dtype=np.uint64)], - ) - csr = read_csr(store, "f", locations=True) - assert csr["locations"].shape == (0,) and csr["locations"].dtype == np.uint64 - - def test_locations_list_length_mismatch_raises(self): - with pytest.raises(ValueError, match="locations_list .* must have the same length"): - write_csr( - MemoryStore(), "f", self._PAYLOADS, self._CELL_IDS, locations_list=self._LOCS[:1] - ) - - def test_per_cell_location_length_mismatch_raises(self): - bad = [self._LOCS[0][:1], self._LOCS[1]] - with pytest.raises(ValueError, match="share the payload's offsets"): - write_csr(MemoryStore(), "f", self._PAYLOADS, self._CELL_IDS, locations_list=bad) - - def test_read_locations_missing_raises(self): - store = MemoryStore() - write_csr(store, "f", self._PAYLOADS, self._CELL_IDS) - with pytest.raises(ValueError, match="has no locations array"): - read_csr(store, "f", locations=True) - - def test_nonempty_locations_on_empty_payload_raise(self): - # An upstream misalignment must fail loudly, not vanish in the filter. - payloads = [np.empty((0, 2), dtype=np.float32)] - locs = [np.array([7], dtype=np.uint64)] - with pytest.raises(ValueError, match="payload is empty"): - write_csr(MemoryStore(), "f", payloads, [0], locations_list=locs) - - def test_read_without_locations_flag_unchanged(self): - store = MemoryStore() - write_csr(store, "f", self._PAYLOADS, self._CELL_IDS, locations_list=self._LOCS) - csr = read_csr(store, "f") - assert set(csr) == {"values", "offsets", "cell_ids"} diff --git a/tests/test_hive.py b/tests/test_hive.py index ad793ee..c1d3f55 100644 --- a/tests/test_hive.py +++ b/tests/test_hive.py @@ -18,7 +18,7 @@ from zagg import hive from zagg.config import default_config, get_data_vars, validate_config from zagg.grids import HealpixGrid -from zagg.grids.morton import morton_decimal, morton_word +from zagg.grids.morton import morton_decimal @pytest.fixture @@ -289,19 +289,37 @@ def _rec(n): class TestProcessAndWriteHive: """Drive ``hive.process_and_write_hive`` with a fake ``process_shard`` that - streams REAL carriers, so the leaf template, dense write, CSR naming, and - stamp ordering are all exercised against real zarr stores.""" + streams REAL carriers, so the leaf template, dense write, ragged vlen + layout (issue #209), and stamp ordering are all exercised against real + zarr stores.""" def _grid(self, cfg): + # Declare the ragged field the streaming fakes emit, so the leaf + # template carries its vlen-bytes array (issue #209). + cfg.aggregation["variables"].setdefault( + "h", + { + "function": "np.sort", + "source": "h_li", + "kind": "ragged", + "inner_shape": [1], + "dtype": "float32", + "fill_value": 0, + }, + ) return HealpixGrid(parent_order=6, child_order=8, layout="fullsphere", config=cfg) def _carrier(self, grid, shard): + from zagg.config import get_agg_fields, get_output_signature + coords = grid.chunk_coords(shard) n = len(coords["cell_ids"]) + agg = get_agg_fields(grid.config) df = pd.DataFrame( { var: np.zeros(n, dtype=np.int32 if var == "count" else np.float32) for var in get_data_vars(grid.config) + if get_output_signature(agg[var])["kind"] != "ragged" } ) for name, vals in coords.items(): @@ -356,11 +374,12 @@ def test_leaf_written_and_stamped(self, monkeypatch, cfg, tmp_path): np.asarray(grp["cell_ids"][:]), np.asarray(grid.chunk_coords(shard)["cell_ids"]), ) - # CSR subgroup named by the shard label (decimal morton string). - label = morton_decimal(shard) - assert morton_word(label) == shard - sub = zarr.open_group(leaf_store, path=f"{grid.group_path}/h/{label}", mode="r") - assert "values" in sub.array_keys() + # The ragged payload sits in the leaf's vlen-bytes array at its cell + # position (issue #209), as ONE data object. + ragged_arr = zarr.open_group(leaf_store, path=grid.group_path, mode="r")["h"] + np.testing.assert_array_equal(np.frombuffer(ragged_arr[0:1][0], " no coverage visible either (issue #200): the tier-0 # payload rides the stamp, so a torn worker never publishes coverage. assert hive.read_coverage(open_store(leaf)) is None - stale = os.path.join(leaf, grid.group_path, "h", morton_decimal(shard)) - assert os.path.exists(stale) # the torn attempt's CSR subgroup - # Plant a sidecar in the debris: the one leaf object zarr does NOT + # The torn attempt's ragged was accumulated, never written (issue #209). + assert not os.path.exists(os.path.join(leaf, grid.group_path, "h", "c")) + stale = os.path.join(leaf, grid.group_path, "stale-debris") + with open(stale, "w") as fh: + fh.write("torn attempt") + # Plant a sidecar in the debris too: the one leaf object zarr does NOT # own must also fall to the wholesale wipe (PR #208 round 2) — this # goes red if the re-template ever drifts to node-by-node rewrites. hive.write_coverage_sidecar(leaf, b"torn-attempt sidecar") @@ -420,7 +444,7 @@ def torn(g, shard_key, urls, **kwargs): assert os.path.exists(sidecar) # Retry (no ragged this time): same leaf, overwritten wholesale — - # the stale subgroup is GONE — and stamped at the end. + # the planted debris is GONE — and stamped at the end. monkeypatch.setattr(processing, "process_shard", self._streaming_fake(grid)) hive.process_and_write_hive( shard, ["s3://bucket/g1.h5"], grid, {}, root, cfg, store_kwargs={} @@ -500,7 +524,9 @@ def wrapped(*a, **k): processing, "write_dataframe_to_zarr", rec("dense", processing.write_dataframe_to_zarr) ) monkeypatch.setattr( - processing, "write_ragged_to_zarr", rec("ragged", processing.write_ragged_to_zarr) + processing, + "write_ragged_leaf_to_zarr", + rec("ragged", processing.write_ragged_leaf_to_zarr), ) monkeypatch.setattr( hive, "write_coverage_sidecar", rec("sidecar", hive.write_coverage_sidecar) diff --git a/tests/test_lambda_handler.py b/tests/test_lambda_handler.py index 1232f6a..6fc3a70 100644 --- a/tests/test_lambda_handler.py +++ b/tests/test_lambda_handler.py @@ -325,10 +325,10 @@ def boom(*a, **k): class TestProcessEventWriteLoop: """Issue #82 phase 7: the handler drives ``process_shard`` with a ``chunk_results`` sink and writes each chunk's dense region (at its own - block_index) plus its ragged (CSR) companion — the same K>1 write loop the - local runner runs. These mock the writers/store to assert the loop's wiring - (sink passed, per-chunk block index used, ragged persisted) without a real - Zarr store.""" + block_index) plus its ragged vlen payloads at the same block (issue #209) — + the same K>1 write loop the local runner runs. These mock the writers/store + to assert the loop's wiring (sink passed, per-chunk block index used, + ragged persisted) without a real Zarr store.""" def _patch(self, handler_mod, monkeypatch, chunks): """Patch process_shard to fill ``chunk_results`` with ``chunks`` and capture @@ -363,13 +363,11 @@ def fake_process_shard(grid, shard_key, granule_urls, **kwargs): def fake_write_dense(carrier, st, *, grid, chunk_idx): cap["dense"].append(chunk_idx) - def fake_write_ragged(ragged, st, *, grid, shard_key): - cap["ragged"].append((shard_key, ragged)) + def fake_write_ragged(ragged, st, *, grid, chunk_idx): + cap["ragged"].append((chunk_idx, ragged)) - # grid stub: a 1-D companion (single-element block index), exposes - # chunk_grid_shape so _block_index_key is exercised on its real path. - # ``chunks_per_shard`` fixes K (issue #91): the handler now derives the - # K==1-vs-K>1 ragged-key choice from the grid, not the chunk count. + # grid stub: a 1-D companion (single-element block index). + # ``chunks_per_shard`` fixes K (issue #91). grid_stub = MagicMock() grid_stub.group_path = "8" grid_stub.chunk_grid_shape = (4,) @@ -377,10 +375,6 @@ def fake_write_ragged(ragged, st, *, grid, shard_key): # Regular (non-sharded) write path — a MagicMock attr is truthy by default, # so pin it off explicitly (issue #108 routes sharded grids elsewhere). grid_stub.sharded = False - # K==1 ragged subgroups are keyed by the grid's shard label (issue #199); - # a deterministic stub (decimal-string shaped) pins that the handler - # routes through the seam rather than stringifying the raw key itself. - grid_stub.shard_label = lambda key: f"-{int(key)}" monkeypatch.setattr(processing, "process_shard", fake_process_shard) monkeypatch.setattr(grids, "from_config", lambda *a, **k: grid_stub) @@ -392,7 +386,7 @@ def fake_write_ragged(ragged, st, *, grid, shard_key): def test_k_gt_1_writes_each_chunk_region(self, handler_mod, monkeypatch): """K=3: the sink loop writes 3 dense regions at distinct block indices, and - each chunk's ragged is keyed by its own _block_index_key (not shard_key).""" + each chunk's ragged lands at the same per-chunk block (issue #209).""" chunks = [ ((0,), pd.DataFrame(), {}), ((1,), pd.DataFrame(), {"h_tdigest": ([], [])}), @@ -405,8 +399,8 @@ def test_k_gt_1_writes_each_chunk_region(self, handler_mod, monkeypatch): assert resp["statusCode"] == 200 # One dense write per chunk, at each chunk's own block_index. assert cap["dense"] == [(0,), (1,), (2,)] - # K>1 -> ragged keyed by _block_index_key(block_index) == 0/1/2, NOT shard_key. - assert [k for k, _r in cap["ragged"]] == [0, 1, 2] + # ... and one ragged write per chunk at the SAME block index. + assert [k for k, _r in cap["ragged"]] == [(0,), (1,), (2,)] def test_non_sharded_streams_via_write_chunk_not_chunk_results(self, handler_mod, monkeypatch): """Issue #91: the non-sharded handler hands ``process_shard`` a ``write_chunk`` @@ -458,10 +452,9 @@ def fake_process_shard(grid, shard_key, granule_urls, **kwargs): assert seen["chunk_results"] is None # no accumulation sink assert written == [(0,), (1,)] # each chunk written as it streamed - def test_k_eq_1_ragged_keyed_by_shard_label(self, handler_mod, monkeypatch): - """K=1: the lone chunk's ragged CSR is persisted, keyed by the grid's - shard label (the decimal morton string for HEALPix — issue #199), not - the raw shard-key int.""" + def test_k_eq_1_ragged_written_at_chunk_block(self, handler_mod, monkeypatch): + """K=1: the lone chunk's ragged payloads are persisted at the chunk's own + block index — the same block the dense write uses (issue #209).""" chunks = [((0,), pd.DataFrame(), {"h_tdigest": ([], [])})] cap = self._patch(handler_mod, monkeypatch, chunks) event = _base_event(_healpix_config_dict()) @@ -469,9 +462,7 @@ def test_k_eq_1_ragged_keyed_by_shard_label(self, handler_mod, monkeypatch): resp = handler_mod._handle_process(event, _context()) assert resp["statusCode"] == 200 assert cap["dense"] == [(0,)] - # Single chunk -> ragged keyed by the stub grid's label of shard_key. - assert len(cap["ragged"]) == 1 - assert cap["ragged"][0][0] == f"-{event['shard_key']}" + assert [k for k, _r in cap["ragged"]] == [(0,)] def test_streaming_later_chunk_write_failure_partial_writes_and_500( self, handler_mod, monkeypatch @@ -842,6 +833,16 @@ class TestProcessHive: def _hive_config_dict(): cfg = default_config("atl06") cfg.output["store_layout"] = "hive" + # A declared ragged field, so the leaf template carries its vlen-bytes + # array (issue #209) and the fakes can stream ragged payloads for it. + cfg.aggregation["variables"]["h"] = { + "function": "np.sort", + "source": "h_li", + "kind": "ragged", + "inner_shape": [1], + "dtype": "float32", + "fill_value": 0, + } return asdict(cfg) def _event(self, tmp_path): @@ -862,14 +863,16 @@ def _grid(): def _carrier(grid, shard): import numpy as np - from zagg.config import get_data_vars + from zagg.config import get_agg_fields, get_data_vars, get_output_signature coords = grid.chunk_coords(shard) n = len(coords["cell_ids"]) + agg = get_agg_fields(grid.config) df = pd.DataFrame( { var: np.zeros(n, dtype=np.int32 if var == "count" else np.float32) for var in get_data_vars(grid.config) + if get_output_signature(agg[var])["kind"] != "ragged" } ) for name, vals in coords.items(): @@ -918,14 +921,19 @@ def test_hive_event_writes_leaf_data_and_stamp(self, handler_mod, monkeypatch, t assert hive.shard_leaf_path(event["store_path"], self._WORD) == expected leaf_store = open_store(expected) - # Dense data landed at the leaf-local block; CSR is label-named. + # Dense data landed at the leaf-local block; the ragged payload sits in + # its vlen-bytes array at the cell position (issue #209), ONE object. grp = zarr.open_group(leaf_store, path=grid.group_path, mode="r", zarr_format=3) np.testing.assert_array_equal( np.asarray(grp["cell_ids"][:]), np.asarray(grid.chunk_coords(self._WORD)["cell_ids"]), ) - sub = zarr.open_group(leaf_store, path=f"{grid.group_path}/h/-4211322", mode="r") - assert "values" in sub.array_keys() + ragged_arr = zarr.open_array(leaf_store, path=f"{grid.group_path}/h", mode="r") + np.testing.assert_array_equal(np.frombuffer(ragged_arr[0:1][0], "1 object. + prefix = f"{self.CHILD}/h_ragged/c/" + assert len([k for k in s_default._store_dict if k.startswith(prefix)]) == 1 + assert len([k for k in s_split._store_dict if k.startswith(prefix)]) > 1 def test_invalid_shard_order_rejected(self): cfg = default_config() @@ -802,11 +804,12 @@ def test_ragged_scalar_vector_coexist(self): assert result["h_raw"].shape == (3, 1) -class TestRaggedCsrWrite: - """Issue #48 phase 4b: cell-resolution ragged (CSR) fields are threaded out of - ``process_shard`` via ``ragged_out`` and persisted by ``write_ragged_to_zarr``, - then read back through the standard ``read_csr`` layout the tensor reader - consumes (``{group_path}/{field}/{shard_key}/values|offsets|cell_ids``).""" +class TestRaggedVlenWrite: + """Cell-resolution ragged fields are threaded out of ``process_shard`` via + ``ragged_out`` (issue #48) and persisted by ``write_ragged_to_zarr`` into + the field's ONE vlen-bytes array (issue #209 — the sharded vlen-bytes + layout replacing the per-shard CSR subgroups): each populated cell holds + the raw little-endian bytes of its ``(n, *inner_shape)`` payload.""" @staticmethod def _ragged_cfg(): @@ -907,13 +910,13 @@ def test_ragged_out_none_is_byte_identical(self, monkeypatch): assert len(result) == 2 assert isinstance(result[0], pd.DataFrame) - def test_end_to_end_write_then_read_csr(self, monkeypatch): - """Full path: process_shard → write_ragged_to_zarr → read_csr returns the - per-cell payloads at the ``{group_path}/{field}/{shard_key}`` CSR layout.""" + def test_end_to_end_write_then_read_vlen(self, monkeypatch): + """Full path: process_shard → write_ragged_to_zarr → the vlen array at + ``{group_path}/{field}`` holds each populated cell's raw payload bytes + at its cell position; empty cells read the ``b""`` fill.""" + import zarr from zarr.storage import MemoryStore - from zagg.csr import iter_csr_cells, read_csr - cfg = self._ragged_cfg() grid, shard_key = self._shard_grid(cfg) children = grid.children(shard_key) @@ -926,30 +929,32 @@ def test_end_to_end_write_then_read_csr(self, monkeypatch): ) self._patch_reads(monkeypatch, df) - # Emit the real product template first: a ragged field must NOT get a - # dense array at ``{group_path}/h_raw`` (which would make the CSR per-shard - # child groups collide with an array node). The dense scalar h_min still - # gets its array. + # The template carries the ragged field as ONE vlen array next to the + # dense scalar (issue #209 — no CSR group prefix anymore), with the + # self-describing element attrs. store = MemoryStore() grid.emit_template(store) - import zarr product = zarr.open_group(store, path=grid.group_path, mode="r") assert "h_min" in product.array_keys() - assert "h_raw" not in product.array_keys() + assert "h_raw" in product.array_keys() + assert product["h_raw"].attrs["ragged"] == { + "spec": "zagg-ragged/1", + "element": {"dtype": "float32", "shape": [-1, 1]}, + } ragged: dict = {} process_shard(grid, shard_key, ["s3://x"], s3_credentials={}, config=cfg, ragged_out=ragged) - write_ragged_to_zarr(ragged, store, grid=grid, shard_key=shard_key) + write_ragged_to_zarr(ragged, store, grid=grid, chunk_idx=grid.block_index(shard_key)) - csr = read_csr(store, f"{grid.group_path}/h_raw/{shard_key}") - cells = dict(iter_csr_cells(csr)) - # Cell positions 1 and 3 are populated; their payloads are the sorted h_li. - assert sorted(cells) == [1, 3] - np.testing.assert_array_equal(cells[1].reshape(-1), [4.0, 5.0]) - np.testing.assert_array_equal(cells[3].reshape(-1), [7.0]) - # The values array carries the declared dtype. - assert csr["values"].dtype == np.dtype("float32") + arr = zarr.open_array(store, path=f"{grid.group_path}/h_raw", mode="r") + base = grid.block_index(shard_key)[0] * grid.cells_per_chunk + block = arr[base : base + grid.cells_per_chunk] + # Cell positions 1 and 3 are populated; their payloads are the sorted + # h_li as raw little-endian float32 bytes (the declared dtype). + assert [i for i, b in enumerate(block) if len(b)] == [1, 3] + np.testing.assert_array_equal(np.frombuffer(block[1], " 0 for w in members: @@ -1115,7 +1136,7 @@ def test_end_to_end_located_write_then_read(self, monkeypatch): assert covered == len(leaf) # every observation is located somewhere # delta=1 forced real merges: at least one location is a coarsened # (below-order-29) enclosing cell, not a raw point word. - all_locs = np.concatenate([locs for _, _, locs in out]) + all_locs = np.concatenate([locs for _, locs in out]) assert not set(int(w) for w in all_locs) <= set(int(w) for w in leaf) def test_empty_cell_returns_located_pair(self): @@ -1155,7 +1176,7 @@ def test_chunk_resolution_located_write_raises(self): {"h_tdigest": ([payload], [0], [locs])}, MemoryStore(), grid=grid, - shard_key=0, + chunk_idx=(0,), ) def test_unlocated_field_keeps_two_tuple(self, monkeypatch): @@ -1177,244 +1198,342 @@ def test_unlocated_field_keeps_two_tuple(self, monkeypatch): assert len(ragged["h_tdigest"]) == 2 -class TestRaggedWriteFanout: - """Issue #142: the sharded path fans out its K ragged (CSR) subgroup writes over - a bounded thread pool instead of a serial loop. The subgroups target disjoint - prefixes, so the on-disk result is identical to serial -- only the write - scheduling changes. These exercise ``_write_ragged_fanout`` directly (hand-built - payloads, no template emit) plus one end-to-end pass through - ``write_shard_to_zarr``.""" - - @staticmethod - def _grid_and_writes(n): - # A real ragged grid supplies group_path + config for write_ragged_to_zarr; - # the payloads are hand-built (distinct key + value per subgroup) so no - # emit_template / process_shard is needed and the test stays fast. - cfg = TestRaggedCsrWrite._ragged_cfg() - grid = HealpixGrid(4, 6, layout="fullsphere", config=cfg) - writes = [ - ({"h_raw": ([np.array([float(k)], dtype=np.float32)], [0])}, k) for k in range(1, n + 1) - ] - return grid, writes +class TestRaggedVlenLayout: + """Issue #209 gates for the sharded vlen-bytes ragged layout (which deleted + the per-inner-chunk CSR fan-out and its issue #142 thread pool / issue #186 + failure roll-call): one object per shard, bit-exact round trip through a + fresh reopen, empty-inner-chunk omission in the shard index, the golden + wire-format pin, the 2-GET single-cell read, and flat-sharded vs hive-leaf + payload parity.""" - def test_fanout_writes_all_subgroups(self): - """Every subgroup lands, even past the 128-worker cap (real concurrent writes).""" - from zagg.csr import read_csr - from zagg.processing.write import _write_ragged_fanout + #: ShardingCodec empty-inner-chunk index marker (zarr v3 spec). + SENTINEL = 2**64 - 1 - grid, writes = self._grid_and_writes(200) # > _RAGGED_WRITE_CONCURRENCY - store = MemoryStore() - _write_ragged_fanout([(r, k) for r, k in writes], store, grid=grid) - for _ragged, key in writes: - csr = read_csr(store, f"{grid.group_path}/h_raw/{key}") - np.testing.assert_array_equal(csr["values"].reshape(-1), [float(key)]) - - def test_fanout_byte_identical_parallel_vs_serial(self, monkeypatch): - """Concurrent fan-out yields byte-for-byte the same store as the serial loop.""" - import zagg.processing.write as wmod - from zagg.processing.write import _write_ragged_fanout - - grid, writes = self._grid_and_writes(6) - # Pin the fd ceiling high so the parallel run truly takes the pool branch - # (else a low-ulimit host could clamp workers to 1 and silently compare - # serial-vs-serial, proving nothing). - monkeypatch.setattr(wmod, "fd_safe_max_workers", lambda: 128) - s_par = MemoryStore() - _write_ragged_fanout([(r, k) for r, k in writes], s_par, grid=grid) - - monkeypatch.setattr(wmod, "_RAGGED_WRITE_CONCURRENCY", 1) # force serial branch - s_ser = MemoryStore() - _write_ragged_fanout([(r, k) for r, k in writes], s_ser, grid=grid) - - assert set(s_par._store_dict) == set(s_ser._store_dict) - for key, val in s_par._store_dict.items(): - assert val.to_bytes() == s_ser._store_dict[key].to_bytes(), f"differ at {key}" - - def test_fanout_cap_respects_fd_ceiling(self, monkeypatch): - """The pool is sized ``min(128, len(writes), fd_safe_max_workers())``; an fd - ceiling of 1 forces the serial branch (no pool constructed).""" - import zagg.processing.write as wmod - from zagg.processing.write import _write_ragged_fanout - - grid, writes = self._grid_and_writes(10) - recorded: dict = {} - real_tpe = wmod.ThreadPoolExecutor - - def spy_tpe(max_workers): - recorded["workers"] = max_workers - return real_tpe(max_workers=max_workers) - - monkeypatch.setattr(wmod, "ThreadPoolExecutor", spy_tpe) - - # fd ceiling (4) clamps below both 128 and len(writes)=10. - monkeypatch.setattr(wmod, "fd_safe_max_workers", lambda: 4) - _write_ragged_fanout([(r, k) for r, k in writes], MemoryStore(), grid=grid) - assert recorded["workers"] == 4 - - # fd ceiling of 1 -> serial branch, ThreadPoolExecutor never constructed. - recorded.clear() - monkeypatch.setattr(wmod, "fd_safe_max_workers", lambda: 1) - _write_ragged_fanout([(r, k) for r, k in writes], MemoryStore(), grid=grid) - assert "workers" not in recorded - - def test_fanout_surfaces_write_failure(self, monkeypatch, caplog): - """A failure in any subgroup write is re-raised (not silently swallowed). - Issue #186 (ask 1): the chained cause never reached the operator — the - worker envelope and CloudWatch only record the summary string — so the - summary now carries the cause's type + message, and the first failure - is logged with its full traceback.""" - import logging - import traceback - - import zagg.processing.write as wmod - from zagg.processing.write import _write_ragged_fanout - - grid, writes = self._grid_and_writes(5) - # Pin the fd ceiling high so the pool branch runs (the serial branch - # re-raises the bare OSError, not the RuntimeError summary). - monkeypatch.setattr(wmod, "fd_safe_max_workers", lambda: 128) - - def boom(ragged, store, *, grid, shard_key): - if shard_key == 3: - raise OSError("disk full") - - monkeypatch.setattr(wmod, "write_ragged_to_zarr", boom) - with caplog.at_level(logging.ERROR, logger="zagg.processing.write"): - with pytest.raises( - RuntimeError, match=r"first failing shard_key 3: OSError: disk full" - ): - _write_ragged_fanout([(r, k) for r, k in writes], MemoryStore(), grid=grid) - rec = next(r for r in caplog.records if "subgroup write failed" in r.message) - assert "shard_key 3" in rec.message - assert "OSError: disk full" in rec.message - # Full traceback attached: the log shows WHERE the cause raised, not - # just its repr — the actual S3/zarr frame in a real reproduction. - assert rec.exc_info is not None - tb_text = "".join(traceback.format_exception(*rec.exc_info)) - assert "OSError: disk full" in tb_text - assert "in boom" in tb_text - # A single failure logs no roll call — it would just repeat the header. - assert not any(r.message.startswith(" ragged (CSR) subgroup ") for r in caplog.records) - - def test_fanout_logs_each_failure_bounded(self, monkeypatch, caplog): - """Multiple subgroup failures each get a one-line ``key: Type: msg`` - summary (issue #186: the race signature ACROSS subgroups), capped at - 20 lines, while exactly one record — the first completed failure — - carries the full traceback.""" - import logging - - import zagg.processing.write as wmod - from zagg.processing.write import _write_ragged_fanout - - grid, writes = self._grid_and_writes(25) - # Pin the fd ceiling high so the run truly takes the pool branch (the - # serial branch raises the first bare exception and logs no roll call). - monkeypatch.setattr(wmod, "fd_safe_max_workers", lambda: 128) - - def boom(ragged, store, *, grid, shard_key): - raise OSError(f"reset {shard_key}") - - monkeypatch.setattr(wmod, "write_ragged_to_zarr", boom) - with caplog.at_level(logging.ERROR, logger="zagg.processing.write"): - with pytest.raises(RuntimeError, match=r"failed for 25 of 25"): - _write_ragged_fanout([(r, k) for r, k in writes], MemoryStore(), grid=grid) - lines = [r.message for r in caplog.records] - # Header first (traceback record), then the roll call as its continuation. - assert lines[0].startswith("ragged (CSR) subgroup write failed") - roll_call = [m for m in lines if m.startswith(" ragged (CSR) subgroup ")] - assert len(roll_call) == wmod._RAGGED_FAILURE_LOG_CAP # bounded roll call - assert any("and 5 more subgroup write failure(s)" in m for m in lines) - # Each roll-call line names its own subgroup key and its own error. - for m in roll_call: - key = m.split("subgroup ")[1].split(":")[0] - assert f"reset {key}" in m - assert sum(1 for r in caplog.records if r.exc_info) == 1 - - def test_fanout_empty_is_noop(self): - """No ragged writes -> nothing written, no pool spun up.""" - from zagg.processing.write import _write_ragged_fanout - - cfg = TestRaggedCsrWrite._ragged_cfg() - grid = HealpixGrid(4, 6, layout="fullsphere", config=cfg) - store = MemoryStore() - _write_ragged_fanout([], store, grid=grid) - assert dict(store._store_dict) == {} - - def test_write_shard_to_zarr_routes_through_fanout(self, monkeypatch): - """``write_shard_to_zarr`` routes its ragged writes through the fan-out - exactly once per shard (wiring), keyed per inner chunk. The fan-out's own - behavior with real payloads is covered by the direct tests above; here we - confirm the sharded write invokes it (and only it) for the CSR side, with - one collected write per populated inner chunk keyed by ``_block_index_key``.""" + @staticmethod + def _sharded_setup(seed=7): + """A sharded K=16 grid + one shard's hand-built chunk_results with + variable-length payloads: some cells empty, every 4th inner chunk + entirely empty. Returns the expected ``{(local_chunk, cell): payload}`` + map alongside.""" from mortie import geo2mort - import zagg.processing.write as wmod + cfg = TestRaggedVlenWrite._ragged_cfg() + grid = HealpixGrid(4, 8, layout="fullsphere", config=cfg, chunk_inner=6, sharded=True) + shard_key = int(geo2mort(-78.5, -132.0, order=4)[0]) + shard_block = grid.block_index(shard_key)[0] + rng = np.random.default_rng(seed) + chunk_results: list = [] + payloads: dict = {} + for block, _children in grid.iter_chunks(shard_key): + local = int(block[0]) - shard_block * grid.chunks_per_shard + if local % 4 == 3: + chunk_results.append((block, pd.DataFrame(), {})) # empty inner chunk + continue + cells = sorted(int(c) for c in rng.choice(grid.cells_per_chunk, 3, replace=False)) + # Variable lengths, INCLUDING an empty payload (n=0) per chunk. + vals = [ + rng.standard_normal((n, 1)).astype(np.float32) + for n in (0, int(rng.integers(1, 6)), int(rng.integers(1, 6))) + ] + chunk_results.append((block, pd.DataFrame(), {"h_raw": (vals, cells)})) + for c, v in zip(cells, vals): + if v.size: + payloads[(local, c)] = v + return grid, shard_key, chunk_results, payloads + + def _write(self, grid, shard_key, chunk_results, store=None): from zagg.processing import write_shard_to_zarr - from zagg.processing.write import _block_index_key - # default_config emits a complete dense template (morton/count/... arrays), - # so the sharded dense write succeeds; it has no ragged field, so the - # fan-out is still invoked once with the (empty) collected list -- proving - # the routing without a bespoke ragged product template. - cfg = default_config() - grid = HealpixGrid(4, 6, layout="fullsphere", config=cfg, chunk_inner=5, sharded=True) - shard_key = int(geo2mort(-78.5, -132.0, order=4)[0]) - children = grid.children(shard_key) - c_first, c_last = int(children[0]), int(children[-1]) - df = pd.DataFrame( - { - "h_li": np.array([3.0, 1.0, 7.0], dtype=np.float32), - "s_li": np.array([0.1, 0.1, 0.1], dtype=np.float32), - "leaf_id": np.array([c_first, c_first, c_last], dtype=np.uint64), - } - ) - calls = {"n": 0} + store = store if store is not None else MemoryStore() + grid.emit_template(store) + write_shard_to_zarr(chunk_results, store, grid=grid, shard_key=shard_key) + return store - def one_shot(*a, **k): - calls["n"] += 1 - return df if calls["n"] == 1 else None + def _data_keys(self, store, grid, field="h_raw"): + prefix = f"{grid.group_path}/{field}/" + return [ + k for k in store._store_dict if k.startswith(prefix) and not k.endswith("zarr.json") + ] - monkeypatch.setattr("zagg.processing._read_group", one_shot) - monkeypatch.setattr("zagg.processing.h5coro.H5Coro", lambda *a, **k: object()) - monkeypatch.setattr("zagg.processing._make_url_rewriter", lambda driver: lambda u: u) - # These tests pin worker/aggregation logic, not the read backend: pin - # the hierarchical delegation seam the one_shot stub intercepts (the - # issue #170 default otherwise resolves to inline, which reads through - # its own chunk-aligned path and hits the object() h5obj stub). - from zagg.index.hierarchical import HierarchicalIndex + def test_one_object_per_shard_and_bitexact_roundtrip(self): + """THE issue #209 gate: the whole shard's ragged payloads land in ONE + store object (was ~7 per populated inner chunk as CSR subgroups), and a + FRESH reopen returns every populated cell bit-exact — empty cells (and + the n=0 payloads) read the ``b""`` fill.""" + import zarr - monkeypatch.setattr( - "zagg.processing.worker.index_from_config", lambda cfg: HierarchicalIndex() - ) + grid, shard_key, chunk_results, payloads = self._sharded_setup() + store = self._write(grid, shard_key, chunk_results) + + assert len(self._data_keys(store, grid)) == 1 + + arr = zarr.open_array(store, path=f"{grid.group_path}/h_raw", mode="r") + base = grid.block_index(shard_key)[0] * grid.cells_per_shard + block = arr[base : base + grid.cells_per_shard] + seen = {} + for i in range(grid.chunks_per_shard): + for c in range(grid.cells_per_chunk): + raw = block[i * grid.cells_per_chunk + c] + if len(raw): + seen[(i, c)] = np.frombuffer(raw, "`` dtype (byte-compatible with numcodecs + ``VLenArray``) without rewriting data.""" + import struct + + from numcodecs import Zstd + + grid, shard_key, _chunk_results, _payloads = self._sharded_setup() + shard_block = grid.block_index(shard_key)[0] + blocks = [b for b, _c in grid.iter_chunks(shard_key)] + first = next(b for b in blocks if int(b[0]) == shard_block * grid.chunks_per_shard) + + p1 = np.array([[1.0], [2.0]], dtype=" 8 bytes + p3 = np.array([[3.0]], dtype=" 4 bytes + chunk_results = [(first, pd.DataFrame(), {"h_raw": ([p1, p3], [1, 3])})] + store = self._write(grid, shard_key, chunk_results) + + (key,) = self._data_keys(store, grid) + shard_bytes = store._store_dict[key].to_bytes() + k_inner = grid.chunks_per_shard + index = np.frombuffer(shard_bytes[-(16 * k_inner + 4) : -4], dtype=" 0) == (len(raw_loc) > 0) # channels populate together + if len(raw): + v = np.frombuffer(raw, "1 with a resolution: chunk scalar companion AND a cell-resolution ragged - field: each chunk writes its own companion slice + CSR group.""" + field: each chunk writes its own companion slice + its own region of the + ragged vlen array (issue #209 — one object per inner chunk unsharded).""" import zarr from mortie import geo2mort from zarr.storage import MemoryStore from zagg.config import default_config - from zagg.csr import read_csr from zagg.grids import from_config from zagg.processing import write_dataframe_to_zarr, write_ragged_to_zarr @@ -1811,21 +1937,29 @@ def test_k_gt_1_with_chunk_companion_and_ragged(self, monkeypatch): grid, shard_key, ["s3://x"], s3_credentials={}, config=cfg, chunk_results=results ) assert len(results) == 4 - n_csr_groups = 0 + n_ragged_chunks = 0 + h_raw = zarr.open_array(store, path=f"{grid.group_path}/h_raw", mode="r") for block_index, carrier, ragged in results: write_dataframe_to_zarr(carrier, store, grid=grid, chunk_idx=block_index) - # Each chunk's ragged keyed by its own block index (K>1). - key = int(block_index[0]) - write_ragged_to_zarr(ragged, store, grid=grid, shard_key=key) + # Each chunk's ragged lands at its own block of the vlen array (K>1). + write_ragged_to_zarr(ragged, store, grid=grid, chunk_idx=block_index) if ragged.get("h_raw") and ragged["h_raw"][0]: - csr = read_csr(store, f"{grid.group_path}/h_raw/{key}") - assert csr["values"].size > 0 - n_csr_groups += 1 + base = int(block_index[0]) * grid.cells_per_chunk + block = h_raw[base : base + grid.cells_per_chunk] + assert sum(len(b) for b in block) > 0 + n_ragged_chunks += 1 # offset_h companion: 4 distinct chunk slices populated. offset = zarr.open_array(store, path="8/offset_h", mode="r")[:] assert int(np.count_nonzero(~np.isnan(offset))) == 4 - # Each chunk had one populated cell -> one CSR group per chunk. - assert n_csr_groups == 4 + # Each chunk had one populated cell -> each chunk's region populated, + # one object per inner chunk on this unsharded path. + assert n_ragged_chunks == 4 + ragged_objs = [ + k + for k in store._store_dict + if k.startswith(f"{grid.group_path}/h_raw/") and not k.endswith("zarr.json") + ] + assert len(ragged_objs) == 4 def test_chunk_precompute_is_per_chunk_not_shard_pooled(self, monkeypatch): """Issue #82 phase 6: a ``chunk_precompute`` anchor is reduced over EACH @@ -2184,7 +2318,6 @@ def test_scalar_vector_ragged_chunk_companions_roundtrip(self, monkeypatch): from zarr.storage import MemoryStore from zagg.config import default_config - from zagg.csr import iter_csr_cells, read_csr from zagg.grids import from_config from zagg.processing import write_dataframe_to_zarr, write_ragged_to_zarr @@ -2244,7 +2377,7 @@ def test_scalar_vector_ragged_chunk_companions_roundtrip(self, monkeypatch): assert len(results) == 1 # K == 1 (no chunk_inner) block_index, carrier, ragged = results[0] write_dataframe_to_zarr(carrier, store, grid=grid, chunk_idx=block_index) - write_ragged_to_zarr(ragged, store, grid=grid, shard_key=shard_key) + write_ragged_to_zarr(ragged, store, grid=grid, chunk_idx=block_index) chunk_idx = grid.block_index(shard_key) # scalar companion: one value at this chunk == min(h_li) == 3.0. @@ -2253,10 +2386,10 @@ def test_scalar_vector_ragged_chunk_companions_roundtrip(self, monkeypatch): # vector companion: the 3-vector edges at this chunk. edges = zarr.open_array(store, path="8/edges_h", mode="r") np.testing.assert_array_equal(edges[chunk_idx], [0.0, 5.0, 10.0]) - # ragged companion: one CSR payload per chunk == edges. - cells = dict(iter_csr_cells(read_csr(store, f"8/edges_ragged/{shard_key}"))) - assert list(cells) == [0] - np.testing.assert_array_equal(cells[0].reshape(-1), [0.0, 5.0, 10.0]) + # ragged companion: one vlen payload at this chunk's block == edges. + ragged_arr = zarr.open_array(store, path="8/edges_ragged", mode="r") + (raw,) = ragged_arr[chunk_idx[0] : chunk_idx[0] + 1] + np.testing.assert_array_equal(np.frombuffer(raw, "1). +""" import math -import warnings import numpy as np +import pandas as pd import pytest import zarr from zarr.storage import MemoryStore -from zagg.csr import write_csr +from zagg.config import PipelineConfig +from zagg.grids import HealpixGrid from zagg.grids.morton import morton_word +from zagg.processing import write_ragged_to_zarr, write_shard_to_zarr from zagg.readers.tdigest_tensor import ( - _resolve_chunk_morton, chunk_z_range, rasterize_cell, + read_cell, + read_locations, read_raw_values, read_tensors, ) from zagg.stats.tdigest import build_tdigest +# Two order-6 shards (decimal morton ids). At K==1 the read chunk IS the +# shard, so the readers report these packed words as the chunk morton ids. +_KEY_A, _KEY_B = "1121121", "2431123" + + +def _cfg(located=False): + """Minimal config: one ragged t-digest field + the morton coordinate.""" + meta = { + "function": "zagg.stats.tdigest.build_tdigest", + "source": "h", + "kind": "ragged", + "inner_shape": [2], + "dtype": "float32", + "fill_value": 0, + } + if located: + meta["location"] = "leaf_id" + return PipelineConfig( + data_source={"groups": ["g"]}, + aggregation={ + "coordinates": {"morton": {"dtype": "uint64", "fill_value": 0}}, + "variables": {"h_tdigest": meta}, + }, + output={"grid": {"type": "healpix", "parent_order": 6, "child_order": 12}}, + ) + + +def _grid(located=False): + """Fullsphere K==1 grid: one 4096-cell (64×64) read chunk per shard.""" + return HealpixGrid(6, 12, layout="fullsphere", config=_cfg(located)) -def _write_chunk(store, field, morton_key, cell_to_values, *, delta=512): - """Write one shard subgroup of per-cell t-digests under {field}/{morton_key}. - ``morton_key`` is the subgroup name: a decimal morton string for the - shard-keyed layout (issue #199), or a plain int for the K>1 block-keyed - layout. - """ +def _write_shard(grid, store, morton_key, cell_to_values, *, delta=512, locations=None): + """Write one shard's digests (and morton coordinate) via the production + per-chunk writer. ``locations`` maps cell -> per-obs uint64 words for a + located field.""" + word = morton_word(morton_key) if isinstance(morton_key, str) else int(morton_key) + block = grid.block_index(word) + base = block[0] * grid.cells_per_chunk + morton_arr = zarr.open_array(store, path=f"{grid.group_path}/morton", mode="r+") + morton_arr[base : base + grid.cells_per_chunk] = grid.children(word) cell_ids = sorted(cell_to_values) - payloads = [build_tdigest(np.asarray(cell_to_values[c]), delta=delta) for c in cell_ids] - write_csr(store, f"{field}/{morton_key}", payloads, cell_ids) + if locations is None: + payloads = [build_tdigest(np.asarray(cell_to_values[c]), delta=delta) for c in cell_ids] + ragged = {"h_tdigest": (payloads, cell_ids)} + else: + pairs = [ + build_tdigest(np.asarray(cell_to_values[c]), delta=delta, locations=locations[c]) + for c in cell_ids + ] + ragged = {"h_tdigest": ([p[0] for p in pairs], cell_ids, [p[1] for p in pairs])} + write_ragged_to_zarr(ragged, store, grid=grid, chunk_idx=block) + return word + + +def _build_store(shards, *, delta=512, located_locs=None): + """Template + shards on the K==1 layout; returns ``(store, grid, words)``.""" + grid = _grid(located=located_locs is not None) + store = MemoryStore() + grid.emit_template(store) + words = {} + for key, cell_to_values in shards.items(): + locs = located_locs[key] if located_locs is not None else None + words[key] = _write_shard(grid, store, key, cell_to_values, delta=delta, locations=locs) + return store, grid, words + + +class _CountingStore(MemoryStore): + """MemoryStore recording every satisfied GET as ``(key, byte_range, nbytes)``.""" + + def __init__(self, *a, **k): + super().__init__(*a, **k) + self.gets: list = [] + + def with_read_only(self, read_only=True): + s = _CountingStore(store_dict=self._store_dict, read_only=read_only) + s.gets = self.gets + return s + + async def get(self, key, prototype, byte_range=None): + r = await super().get(key, prototype, byte_range) + if r is not None: + self.gets.append((key, byte_range, len(r))) + return r + + +def _sharded_store(store, *, populate): + """A K=16 ShardingCodec'd store with one populated cell per chunk in + ``populate`` (local chunk ordinals); returns ``(grid, shard_word, + {local: global_cell})``.""" + rng = np.random.default_rng(21) + cfg = _cfg() + cfg.output["grid"]["chunk_inner"] = 8 + cfg.output["grid"]["sharded"] = True + grid = HealpixGrid(6, 12, layout="fullsphere", config=cfg, chunk_inner=8, sharded=True) + shard6 = morton_word(_KEY_A) + grid.emit_template(store) + morton_arr = zarr.open_array(store, path="12/morton", mode="r+") + shard_block = grid.block_index(shard6)[0] + chunk_results, targets = [], {} + for block, children in grid.iter_chunks(shard6): + local = int(block[0]) - shard_block * grid.chunks_per_shard + base = int(block[0]) * grid.cells_per_chunk + morton_arr[base : base + grid.cells_per_chunk] = np.asarray(children) + if local not in populate: + chunk_results.append((block, pd.DataFrame(), {})) + continue + payloads = [build_tdigest(rng.uniform(0.0, 30.0, 300), delta=512)] + chunk_results.append((block, pd.DataFrame(), {"h_tdigest": (payloads, [11])})) + targets[local] = base + 11 + write_shard_to_zarr(chunk_results, store, grid=grid, shard_key=shard6) + return grid, shard6, targets class TestRasterizeCell: @@ -170,47 +282,36 @@ def test_unknown_fit_raises(self): class TestReadTensors: - # Two coverage chunks named by their decimal morton strings (issue #199); - # the readers report the parsed packed words. - _KEY_A, _KEY_B = "112", "243" - def _store(self): - store = MemoryStore() rng = np.random.default_rng(10) - # Two chunks (decimal morton names), a few populated cells each. - _write_chunk( - store, - "h_tdigest", - self._KEY_A, + store, _grid_, _words = _build_store( { - 0: rng.uniform(10.0, 30.0, 3_000), - 5: rng.uniform(12.0, 28.0, 2_000), - 4095: rng.uniform(11.0, 29.0, 1_500), - }, - ) - _write_chunk( - store, - "h_tdigest", - self._KEY_B, - {7: rng.uniform(40.0, 60.0, 2_500), 63: rng.uniform(42.0, 58.0, 2_000)}, + _KEY_A: { + 0: rng.uniform(10.0, 30.0, 3_000), + 5: rng.uniform(12.0, 28.0, 2_000), + 4095: rng.uniform(11.0, 29.0, 1_500), + }, + _KEY_B: {7: rng.uniform(40.0, 60.0, 2_500), 63: rng.uniform(42.0, 58.0, 2_000)}, + } ) return store def test_shape_and_dtype_default(self): - out = dict((m, t) for t, m in read_tensors(self._store(), "h_tdigest")) - assert set(out) == {morton_word(self._KEY_A), morton_word(self._KEY_B)} + out = dict((m, t) for t, m in read_tensors(self._store(), "12/h_tdigest")) + assert set(out) == {morton_word(_KEY_A), morton_word(_KEY_B)} for t in out.values(): assert t.shape == (64, 64, 128) assert t.dtype == np.uint32 - def test_morton_recovered_from_subgroup_name(self): - # Round trip at the read boundary: decimal subgroup name -> packed word. - mortons = sorted(m for _, m in read_tensors(self._store(), "h_tdigest")) - assert mortons == sorted(morton_word(k) for k in (self._KEY_A, self._KEY_B)) + def test_morton_derived_from_coordinate(self): + # Chunk identity round trip: the per-cell morton coordinate coarsens + # to the chunk's coverage-cell id (the packed shard word at K==1). + mortons = sorted(m for _, m in read_tensors(self._store(), "12/h_tdigest")) + assert mortons == sorted(morton_word(k) for k in (_KEY_A, _KEY_B)) def test_populated_cell_placement_rowmajor(self): - out = dict((m, t) for t, m in read_tensors(self._store(), "h_tdigest")) - t = out[morton_word(self._KEY_A)] + out = dict((m, t) for t, m in read_tensors(self._store(), "12/h_tdigest")) + t = out[morton_word(_KEY_A)] # cell 5 → row 0, col 5; cell 4095 → row 63, col 63. assert t[0, 5].sum() > 0 assert t[63, 63].sum() > 0 @@ -218,12 +319,11 @@ def test_populated_cell_placement_rowmajor(self): assert t[10, 10].sum() == 0 def test_counts_match_population(self): - store = MemoryStore() rng = np.random.default_rng(11) n = 5_000 - _write_chunk(store, "f", "11", {0: rng.uniform(0.0, 40.0, n)}) - t, m = next(read_tensors(store, "f", n_bins=128, resolution=0.5)) - assert m == morton_word("11") + store, _g, words = _build_store({_KEY_A: {0: rng.uniform(0.0, 40.0, n)}}) + t, m = next(read_tensors(store, "12/h_tdigest", n_bins=128, resolution=0.5)) + assert m == words[_KEY_A] # Most of the population should land in-window (uniform [0,40] in a 64 m # window anchored at floor of the 5th pct). assert 0.8 * n <= t[0, 0].sum() <= n @@ -233,124 +333,174 @@ def test_counts_match_population(self): [("uint16", np.uint16), ("uint32", np.uint32), ("float32", np.float32)], ) def test_dtype_flag(self, dtype, np_dtype): - store = MemoryStore() rng = np.random.default_rng(12) - _write_chunk(store, "f", 1, {0: rng.uniform(0.0, 30.0, 2_000)}) - t, _ = next(read_tensors(store, "f", dtype=dtype)) + store, _g, _w = _build_store({_KEY_A: {0: rng.uniform(0.0, 30.0, 2_000)}}) + t, _ = next(read_tensors(store, "12/h_tdigest", dtype=dtype)) assert t.dtype == np_dtype - def test_morton_coord_array_preferred(self): - # Block-keyed (K>1) subgroup names — "100"/"250" are not decimal morton - # ids, so the reader falls back to the int parse for the whole store. - store = MemoryStore() - rng = np.random.default_rng(13) - _write_chunk(store, "f", 100, {0: rng.uniform(0.0, 30.0, 2_000)}) - _write_chunk(store, "f", 250, {0: rng.uniform(0.0, 30.0, 2_000)}) - # Sibling coord maps sorted chunk order [100, 250] → custom mortons. - arr = zarr.open_array( - store, path="f/morton", mode="w", shape=(2,), chunks=(2,), dtype="uint64" - ) - arr[...] = np.array([900, 901], dtype=np.uint64) - mortons = sorted(m for _, m in read_tensors(store, "f")) - assert mortons == [900, 901] - - def test_morton_coord_array_mixed_digit_keys(self): - # Regression: subgroup names of differing digit counts must align with - # the coord array in numeric (not lexicographic) order — a lexicographic - # zip would pair "1000" before "99" and mis-assign mortons. - store = MemoryStore() - rng = np.random.default_rng(33) - for key in (99, 100, 1000): - _write_chunk(store, "f", key, {0: rng.uniform(0.0, 30.0, 1_000)}) - arr = zarr.open_array( - store, path="f/morton", mode="w", shape=(3,), chunks=(3,), dtype="uint64" - ) - # Coord in ascending-morton chunk order (99, 100, 1000). - arr[...] = np.array([900, 901, 902], dtype=np.uint64) - mapping = _resolve_chunk_morton(store, "f", ["99", "100", "1000"], 3) - assert mapping == {"99": 900, "100": 901, "1000": 902} - - def test_morton_coord_array_decimal_names_word_order(self): - # Decimal morton names align with the coord array in packed-word order - # (issue #199): northern-base "31123" sorts before southern "-31123", - # even though int/lexicographic order would put the negative first. - store = MemoryStore() - rng = np.random.default_rng(34) - for key in ("-31123", "31123"): - _write_chunk(store, "f", key, {0: rng.uniform(0.0, 30.0, 1_000)}) - assert morton_word("31123") < morton_word("-31123") - arr = zarr.open_array( - store, path="f/morton", mode="w", shape=(2,), chunks=(2,), dtype="uint64" - ) - arr[...] = np.array([900, 901], dtype=np.uint64) - mapping = _resolve_chunk_morton(store, "f", ["-31123", "31123"], 3) - assert mapping == {"31123": 900, "-31123": 901} - - def test_lone_digit_names_are_block_keys(self): - # Review finding (PR #205): "5"/"6" fit the raw decimal grammar as - # order-0 ids, but real shard ids always carry >= 2 digits (shard - # orders are >= 1) — lone digits are block-index keys, parsed as ints. - from zagg.readers.tdigest_tensor import _subgroup_words - - assert _subgroup_words(["5", "6", "1"]) == {"5": 5, "6": 6, "1": 1} - - def test_signed_name_in_block_keyed_store_raises(self): - # Review finding (PR #205): a signed name can never be a block index, - # so a store mixing decimal-id and non-id names (pre-#199 debris) must - # fail loudly instead of silently degrading every id to a bogus int. - from zagg.readers.tdigest_tensor import _subgroup_words - - with pytest.raises(ValueError, match="ambiguous"): - _subgroup_words(["-31123", "100"]) - def test_raise_when_chunk_too_wide(self): - store = MemoryStore() rng = np.random.default_rng(14) - _write_chunk(store, "f", 1, {0: rng.uniform(0.0, 400.0, 5_000)}) + store, _g, _w = _build_store({_KEY_A: {0: rng.uniform(0.0, 400.0, 5_000)}}) with pytest.raises(ValueError, match="exceeds the fixed window"): - next(read_tensors(store, "f", bottom=0.0, top=1.0)) + next(read_tensors(store, "12/h_tdigest", bottom=0.0, top=1.0)) def test_degrade_resolution_fits(self): - store = MemoryStore() rng = np.random.default_rng(15) - _write_chunk(store, "f", 1, {0: rng.uniform(0.0, 400.0, 5_000)}) - t, _ = next(read_tensors(store, "f", bottom=0.0, top=1.0, fit="degrade_resolution")) + store, _g, _w = _build_store({_KEY_A: {0: rng.uniform(0.0, 400.0, 5_000)}}) + t, _ = next( + read_tensors(store, "12/h_tdigest", bottom=0.0, top=1.0, fit="degrade_resolution") + ) assert t.shape == (64, 64, 128) def test_unknown_dtype_raises(self): - store = MemoryStore() rng = np.random.default_rng(16) - _write_chunk(store, "f", 1, {0: rng.uniform(0.0, 30.0, 1_000)}) + store, _g, _w = _build_store({_KEY_A: {0: rng.uniform(0.0, 30.0, 1_000)}}) with pytest.raises(ValueError, match="unknown dtype"): - next(read_tensors(store, "f", dtype="float64")) + next(read_tensors(store, "12/h_tdigest", dtype="float64")) + + def test_missing_element_attrs_is_pointed_error(self): + # Forward-compat guard: a vlen array WITHOUT the writer's element + # declaration must fail loudly, not decode under a guessed layout. + import warnings + + from zarr.codecs import VLenBytesCodec + + store = MemoryStore() + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + a = zarr.create_array( + store, + name="12/plain", + shape=(16,), + chunks=(16,), + dtype="bytes", + serializer=VLenBytesCodec(), + fill_value=b"", + ) + obj = np.empty(16, dtype=object) + obj[:] = b"" + obj[0] = b"\x00\x00\x80?" + a[:] = obj + with pytest.raises(ValueError, match="no ragged element declaration"): + list(read_tensors(store, "12/plain")) + + def test_missing_morton_sibling_is_pointed_error(self): + # The chunk-identity source is the sibling morton coordinate; a store + # without it cannot be swept (hard break, no name fallback). + rng = np.random.default_rng(17) + store, _grid_, _w = _build_store({_KEY_A: {0: rng.uniform(0.0, 30.0, 500)}}) + clone = MemoryStore() + # Copy everything except the morton coordinate array. + clone._store_dict.update( + {k: v for k, v in store._store_dict.items() if not k.startswith("12/morton")} + ) + with pytest.raises(ValueError, match="no sibling 'morton' coordinate"): + list(read_tensors(clone, "12/h_tdigest")) + + def test_unknown_spec_raises_pointed(self): + # Version gate (review, PR #211 round 2): a future writer's attrs must + # be adopted deliberately, never half-parsed — the coverage-envelope + # discipline. Issue #210's typed-dtype migration bumps this seam. + store, _g, _w = _build_store({_KEY_A: {0: np.array([1.0, 2.0])}}) + arr = zarr.open_array(store, path="12/h_tdigest", mode="r+") + arr.attrs["ragged"] = {**arr.attrs["ragged"], "spec": "zagg-ragged/2"} + with pytest.raises(ValueError, match="understands 'zagg-ragged/1' only"): + list(read_tensors(store, "12/h_tdigest")) + + def test_garbage_element_dtype_raises_pointed(self): + # A corrupt dtype string surfaces as this layout's pointed error, not + # a bare numpy TypeError with no mention of the field or contract. + store, _g, _w = _build_store({_KEY_A: {0: np.array([1.0, 2.0])}}) + arr = zarr.open_array(store, path="12/h_tdigest", mode="r+") + meta = dict(arr.attrs["ragged"]) + meta["element"] = {**meta["element"], "dtype": "not-a-dtype"} + arr.attrs["ragged"] = meta + with pytest.raises(ValueError, match="unreadable element dtype"): + list(read_tensors(store, "12/h_tdigest")) + + def test_sharded_and_regular_layouts_read_identically(self): + """Q1 pin (issue #209): the K>1 flat layouts differ on disk — one + object per shard (ShardingCodec) vs one per inner chunk (regular) — + but both are self-describing, so ONE reader path yields identical + tensors + chunk ids from the same logical data.""" + from mortie import generate_morton_children + + rng = np.random.default_rng(18) + shard6 = morton_word(_KEY_A) + + # Grid B: order-6 shards of K=16 order-8 inner chunks, ShardingCodec. + cfg = _cfg() + cfg.output["grid"]["chunk_inner"] = 8 + cfg.output["grid"]["sharded"] = True + grid_b = HealpixGrid(6, 12, layout="fullsphere", config=cfg, chunk_inner=8, sharded=True) + # Grid A: order-8 shards (chunk == shard), regular chunks — the SAME + # 256-cell order-8 read chunks as B's inner chunks. + grid_a = HealpixGrid(8, 12, layout="fullsphere", config=_cfg()) + assert grid_a.cells_per_chunk == grid_b.cells_per_chunk == 256 + + # The same per-cell data, addressed by order-8 sub-shard. + data = { + int(sub): {int(rng.integers(0, 256)): rng.uniform(5.0, 25.0, 400)} + for sub in np.asarray(generate_morton_children(shard6, 8))[[0, 5, 11]] + } + + store_a = MemoryStore() + grid_a.emit_template(store_a) + for sub, cells in data.items(): + _write_shard(grid_a, store_a, sub, cells) + + store_b = MemoryStore() + grid_b.emit_template(store_b) + morton_b = zarr.open_array(store_b, path="12/morton", mode="r+") + # iter_chunks enumerates the shard's order-8 sub-cells in the same + # (ascending morton-children) order generate_morton_children yields. + subs = [int(s) for s in np.asarray(generate_morton_children(shard6, 8))] + chunk_results = [] + for (block, children), sub in zip(grid_b.iter_chunks(shard6), subs): + base = int(block[0]) * grid_b.cells_per_chunk + morton_b[base : base + grid_b.cells_per_chunk] = np.asarray(children) + cells = data.get(sub) + if cells is None: + chunk_results.append((block, pd.DataFrame(), {})) + continue + cell_ids = sorted(cells) + payloads = [build_tdigest(np.asarray(cells[c]), delta=512) for c in cell_ids] + chunk_results.append((block, pd.DataFrame(), {"h_tdigest": (payloads, cell_ids)})) + write_shard_to_zarr(chunk_results, store_b, grid=grid_b, shard_key=shard6) + + out_a = sorted(read_tensors(store_a, "12/h_tdigest"), key=lambda tm: tm[1]) + out_b = sorted(read_tensors(store_b, "12/h_tdigest"), key=lambda tm: tm[1]) + assert [m for _t, m in out_a] == [m for _t, m in out_b] == sorted(data) + for (t_a, _ma), (t_b, _mb) in zip(out_a, out_b): + assert t_a.shape == (16, 16, 128) + np.testing.assert_array_equal(t_a, t_b) class TestReadRawValues: def test_recovers_unmerged_samples_exactly(self): - store = MemoryStore() # Few enough values (< delta) that build_tdigest performs no merges. vals = np.array([3.0, 1.0, 2.0, 5.0, 4.0]) - _write_chunk(store, "f", "42", {7: vals}, delta=512) - out = list(read_raw_values(store, "f")) + store, _g, words = _build_store({_KEY_A: {7: vals}}) + out = list(read_raw_values(store, "12/h_tdigest")) assert len(out) == 1 morton, cell_id, recovered = out[0] - assert morton == morton_word("42") + assert morton == words[_KEY_A] assert cell_id == 7 # Digest stores centroids sorted by mean → sorted samples. np.testing.assert_allclose(recovered, np.sort(vals)) def test_merged_digest_raises(self): - store = MemoryStore() rng = np.random.default_rng(20) # Many values at small delta → merges (weight > 1) somewhere. - _write_chunk(store, "f", 1, {0: rng.standard_normal(5_000)}, delta=64) + store, _g, _w = _build_store({_KEY_A: {0: rng.standard_normal(5_000)}}, delta=64) with pytest.raises(ValueError, match="not losslessly recoverable"): - list(read_raw_values(store, "f")) + list(read_raw_values(store, "12/h_tdigest")) class TestReadLocations: """Issue #87: the location-channel reader yields per-cell uint64 morton - vectors aligned with the digest rows the other readers see.""" + vectors aligned with the digest rows the other readers see, bound through + the payload array's attrs declaration (issue #209).""" @staticmethod def _point_words(n, seed): @@ -359,27 +509,20 @@ def _point_words(n, seed): return point_words(n, seed) def _located_store(self, delta=512): - store = MemoryStore() vals = {3: np.array([3.0, 1.0, 2.0]), 9: np.array([5.0, 4.0])} locs_in = {3: self._point_words(3, 1), 9: self._point_words(2, 2)} - cell_ids = sorted(vals) - pairs = [build_tdigest(vals[c], delta=delta, locations=locs_in[c]) for c in cell_ids] - write_csr( - store, - "f/42", - [p[0] for p in pairs], - cell_ids, - locations_list=[p[1] for p in pairs], - ) - return store, vals, locs_in + store, _g, words = _build_store({_KEY_A: vals}, delta=delta, located_locs={_KEY_A: locs_in}) + return store, vals, locs_in, words[_KEY_A] - def test_yields_per_cell_uint64_vectors(self): - from zagg.readers.tdigest_tensor import read_locations + def test_linkage_declared_in_attrs(self): + store, _vals, _locs, _w = self._located_store() + payload = zarr.open_array(store, path="12/h_tdigest", mode="r") + assert payload.attrs["ragged"]["locations"] == "h_tdigest_locations" - store, vals, locs_in = self._located_store() - out = {(m, c): locs for m, c, locs in read_locations(store, "f")} - w = morton_word("42") - assert set(out) == {(w, 3), (w, 9)} + def test_yields_per_cell_uint64_vectors(self): + store, vals, locs_in, word = self._located_store() + out = {(m, c): locs for m, c, locs in read_locations(store, "12/h_tdigest")} + assert set(out) == {(word, 3), (word, 9)} for (_, cid), locs in out.items(): assert locs.dtype == np.uint64 # Loss-free regime: locations are the cell's point words co-sorted @@ -388,30 +531,91 @@ def test_yields_per_cell_uint64_vectors(self): np.testing.assert_array_equal(locs, expected) def test_aligned_with_read_raw_values(self): - from zagg.readers.tdigest_tensor import read_locations - - store, vals, _ = self._located_store() - raw = {(m, c): v for m, c, v in read_raw_values(store, "f")} - locs = {(m, c): loc for m, c, loc in read_locations(store, "f")} + store, _vals, _locs, _w = self._located_store() + raw = {(m, c): v for m, c, v in read_raw_values(store, "12/h_tdigest")} + locs = {(m, c): loc for m, c, loc in read_locations(store, "12/h_tdigest")} assert set(raw) == set(locs) for key in raw: assert len(raw[key]) == len(locs[key]) def test_value_only_field_raises_clearly(self): - from zagg.readers.tdigest_tensor import read_locations - - store = MemoryStore() - _write_chunk(store, "f", 7, {0: np.array([1.0, 2.0])}) - with pytest.raises(ValueError, match="has no locations array"): - list(read_locations(store, "f")) + store, _g, _w = _build_store({_KEY_A: {0: np.array([1.0, 2.0])}}) + with pytest.raises(ValueError, match="declares no locations channel"): + list(read_locations(store, "12/h_tdigest")) + + +class TestReadCell: + """Issue #209 single-cell path: index the vlen array directly.""" + + def test_decodes_declared_element(self): + vals = np.array([3.0, 1.0, 2.0]) + store, grid, words = _build_store({_KEY_A: {7: vals}}) + base = grid.block_index(words[_KEY_A])[0] * grid.cells_per_chunk + digest = read_cell(store, "12/h_tdigest", base + 7) + assert digest.shape == (3, 2) + np.testing.assert_allclose(digest[:, 0], np.sort(vals)) + # An absent cell decodes to the zero-length element. + assert read_cell(store, "12/h_tdigest", base + 8).shape == (0, 2) + + def test_reads_locations_sibling(self): + vals = {3: np.array([3.0, 1.0, 2.0])} + locs = {3: TestReadLocations._point_words(3, 4)} + store, grid, words = _build_store({_KEY_A: vals}, located_locs={_KEY_A: locs}) + base = grid.block_index(words[_KEY_A])[0] * grid.cells_per_chunk + out = read_cell(store, "12/h_tdigest_locations", base + 3) + assert out.dtype == np.uint64 and out.shape == (3,) + + def test_out_of_range_and_negative_raise_pointed(self): + """Review (PR #211 round 2): zarr basic selection CLAMPS the slice, so + without the bounds guard an out-of-range cell died on a bare unpack; + -1 gets neither numpy wrap nor a comprehensible error.""" + store, _grid_, _w = _build_store({_KEY_A: {0: np.array([1.0, 2.0])}}) + n = zarr.open_array(store, path="12/h_tdigest", mode="r").shape[0] + with pytest.raises(IndexError, match=f"cell {n} out of range"): + read_cell(store, "12/h_tdigest", n) + with pytest.raises(IndexError, match="negative indices do not wrap"): + read_cell(store, "12/h_tdigest", -1) + + def test_two_ranged_gets_on_sharded_store(self): + """One cell from a ShardingCodec'd store = exactly 2 ranged GETs on + the shard object (index suffix + one inner chunk), never the whole + object — the random-access property the layout was chosen for.""" + store = _CountingStore() + grid, _shard6, targets = _sharded_store(store, populate=range(16)) + + store.gets.clear() + digest = read_cell(store, "12/h_tdigest", targets[0]) + assert digest.shape[1] == 2 and digest.shape[0] > 0 + data_gets = [g for g in store.gets if "/h_tdigest/c/" in g[0]] + assert len(data_gets) == 2 + (obj_key, _r0, n0), (_k1, _r1, n1) = data_gets + obj_size = len(store._store_dict[obj_key].to_bytes()) + assert n0 == 16 * grid.chunks_per_shard + 4 # the shard-index suffix + assert n1 <= obj_size // 4 # one compressed inner chunk, not the object + + def test_sweep_reads_each_stored_object_once(self): + """Review (PR #211 round 2): the whole-store sweep reads each stored + span in ONE slice, so a stored object is fetched exactly ONCE (zarr + reads a whole outer chunk in a single GET and splits it locally) — + never a per-inner-chunk index re-fetch (K identical suffix reads per + shard, ~256 redundant round trips at production geometry).""" + store = _CountingStore() + _grid_, _shard6, targets = _sharded_store(store, populate={0, 3, 7, 12}) + assert len(targets) == 4 + + store.gets.clear() + out = list(read_tensors(store, "12/h_tdigest")) + assert len(out) == 4 + data_gets = [g for g in store.gets if "/h_tdigest/c/" in g[0]] + assert len(data_gets) == 1 # the whole stored object, once class TestReadParityWithoutConsolidation: """Issue #191: consolidation is now opt-out, so published stores are read without a consolidated-metadata blob. Pin that every reader navigates a non-consolidated store to the same bytes it would read from a consolidated - one — readers reach paths directly (``zarr.open_group``/``open_array`` per - subgroup), never the consolidated blob.""" + one — readers reach paths directly (``zarr.open_array`` per node), never + the consolidated blob.""" @staticmethod def _point_words(n, seed): @@ -420,52 +624,30 @@ def _point_words(n, seed): return point_words(n, seed) def _build_store(self): - """A two-shard located t-digest field, written like the worker's CSR - path (per-parent-morton subgroups), with no consolidated metadata. - - The leaf ShardingCodec is deliberately omitted: ``consolidate_metadata`` - is a pure metadata-tree op (it gathers each group/array ``zarr.json``), - orthogonal to the leaf codec, and the per-morton CSR subgroups are - exactly the objects that dominate the finalize cost this PR skips — so - the parity claim holds without the sharding codec in the fixture.""" - store = MemoryStore() - # Two coverage chunks (parent mortons 100 and 205), distinct values so - # read_raw_values is lossless (every centroid weight 1). + """A two-shard located t-digest field on the vlen layout (issue #209), + written through the production template + per-chunk writer, with no + consolidated metadata.""" shards = { - 100: {0: np.array([10.0, 11.0, 12.0]), 5: np.array([13.0, 14.0])}, - 205: {2: np.array([20.0, 21.0]), 63: np.array([22.0, 23.0, 24.0])}, + _KEY_A: {0: np.array([10.0, 11.0, 12.0]), 5: np.array([13.0, 14.0])}, + _KEY_B: {2: np.array([20.0, 21.0]), 63: np.array([22.0, 23.0, 24.0])}, } - seed = 1 - for morton, cell_to_vals in shards.items(): - cell_ids = sorted(cell_to_vals) - pairs = [ - build_tdigest( - cell_to_vals[c], - delta=512, - locations=self._point_words(len(cell_to_vals[c]), seed), - ) - for c in cell_ids - ] - write_csr( - store, - f"h_tdigest/{morton}", - [p[0] for p in pairs], - cell_ids, - locations_list=[p[1] for p in pairs], - ) - seed += 1 + locs = { + key: {c: self._point_words(len(v), seed) for c, v in cells.items()} + for seed, (key, cells) in enumerate(shards.items(), start=1) + } + store, _g, _w = _build_store(shards, located_locs=locs) return store @staticmethod def _read_all(store): - from zagg.readers.tdigest_tensor import read_locations - - tensors = {m: t for t, m in read_tensors(store, "h_tdigest", n_bins=64, resolution=0.5)} - raw = {(m, c): v for m, c, v in read_raw_values(store, "h_tdigest")} - locs = {(m, c): loc for m, c, loc in read_locations(store, "h_tdigest")} + tensors = {m: t for t, m in read_tensors(store, "12/h_tdigest", n_bins=64, resolution=0.5)} + raw = {(m, c): v for m, c, v in read_raw_values(store, "12/h_tdigest")} + locs = {(m, c): loc for m, c, loc in read_locations(store, "12/h_tdigest")} return tensors, raw, locs def test_non_consolidated_is_navigable_and_reads_parity(self): + import warnings + store = self._build_store() # A freshly written store carries NO consolidated-metadata blob — the @@ -475,8 +657,9 @@ def test_non_consolidated_is_navigable_and_reads_parity(self): tensors_plain, raw_plain, locs_plain = self._read_all(store) # Sanity: the readers actually reached the data. - assert set(tensors_plain) == {100, 205} - assert (100, 0) in raw_plain and (205, 63) in raw_plain + word_a, word_b = morton_word(_KEY_A), morton_word(_KEY_B) + assert set(tensors_plain) == {word_a, word_b} + assert (word_a, 0) in raw_plain and (word_b, 63) in raw_plain # Consolidate the SAME store and re-read: consolidation only adds a # metadata blob no reader consults, so every byte must match. diff --git a/tests/test_runner.py b/tests/test_runner.py index 254ac0d..e055e2e 100644 --- a/tests/test_runner.py +++ b/tests/test_runner.py @@ -1216,9 +1216,9 @@ class TestProcessAndWriteStreaming: """Issue #91: the non-sharded ``_process_and_write`` streams each chunk through a ``write_chunk`` callback (no ``chunk_results`` accumulation). Drive a fake ``process_shard`` that streams 1 and K>1 chunks through the callback and assert - the dense ``chunk_idx`` sequence + ragged keying (the grid's shard label at K=1 - — issue #199 — block-index key at K>1) — the runner-level analogue of the - lambda streaming test.""" + the dense + ragged writes both land at each chunk's own block index (the + ragged vlen array is block-addressed like the dense arrays — issue #209) — + the runner-level analogue of the lambda streaming test.""" def _run(self, monkeypatch, atl06_config, *, chunks_per_shard, chunks): from unittest.mock import MagicMock @@ -1240,9 +1240,6 @@ def fake_process_shard(grid, shard_key, urls, **kwargs): grid.sharded = False grid.chunks_per_shard = chunks_per_shard grid.chunk_grid_shape = (4,) - # K==1 ragged subgroups are keyed by the grid's shard label (issue - # #199); a deterministic decimal-string-shaped stub pins the seam. - grid.shard_label = lambda key: f"-{int(key)}" monkeypatch.setattr(runner, "process_shard", fake_process_shard) monkeypatch.setattr( @@ -1253,10 +1250,8 @@ def fake_process_shard(grid, shard_key, urls, **kwargs): monkeypatch.setattr( runner, "write_ragged_to_zarr", - lambda r, st, *, grid, shard_key: cap["ragged"].append(shard_key), + lambda r, st, *, grid, chunk_idx: cap["ragged"].append(chunk_idx), ) - # _block_index_key on a 1-D grid is block_index[0]. - monkeypatch.setattr(runner, "_block_index_key", lambda b, g: int(b[0])) runner._process_and_write( 5, (5,), @@ -1269,7 +1264,7 @@ def fake_process_shard(grid, shard_key, urls, **kwargs): ) return cap - def test_k1_streams_ragged_keyed_by_shard_label(self, monkeypatch, atl06_config): + def test_k1_streams_ragged_at_chunk_block(self, monkeypatch, atl06_config): import pandas as pd cap = self._run( @@ -1281,9 +1276,9 @@ def test_k1_streams_ragged_keyed_by_shard_label(self, monkeypatch, atl06_config) # Streaming seam wired: callback passed, no accumulation sink. assert callable(cap["write_chunk"]) and cap["chunk_results"] is None assert cap["dense"] == [(5,)] - assert cap["ragged"] == ["-5"] # K=1 -> keyed by the grid's shard label (#199) + assert cap["ragged"] == [(5,)] # same block as the dense write (issue #209) - def test_k_gt_1_streams_ragged_keyed_by_block_index(self, monkeypatch, atl06_config): + def test_k_gt_1_streams_ragged_at_each_chunk_block(self, monkeypatch, atl06_config): import pandas as pd cap = self._run( @@ -1297,7 +1292,7 @@ def test_k_gt_1_streams_ragged_keyed_by_block_index(self, monkeypatch, atl06_con ], ) assert cap["dense"] == [(0,), (1,), (2,)] - assert cap["ragged"] == [0, 1, 2] # K>1 -> keyed by _block_index_key + assert cap["ragged"] == [(0,), (1,), (2,)] # K>1 -> each chunk's own block def _stub_grid(): @@ -1308,7 +1303,7 @@ def _stub_grid(): grid.spatial_signature.return_value = {} grid.block_index.side_effect = lambda k: (k,) # Deterministic decimal-string-shaped labels (issue #199) so status keys - # and CSR subgroup names built off the stub are real strings. + # built off the stub are real strings. grid.shard_label.side_effect = lambda k: f"-{int(k)}" grid.emit_template.side_effect = lambda store, overwrite=False: store return grid