Skip to content

notebooks: binder-runnable reader demos (hhdc_viewer, waveform_viewer) - #540

Merged
espg merged 44 commits into
mainfrom
claude/binder-reader-notebooks
Aug 29, 2026
Merged

notebooks: binder-runnable reader demos (hhdc_viewer, waveform_viewer)#540
espg merged 44 commits into
mainfrom
claude/binder-reader-notebooks

Conversation

@espg

@espg espg commented Aug 26, 2026

Copy link
Copy Markdown
Member

Replaces #328 (closed 2026-08-27).

Merge as a SQUASH. Two .npy blobs (6 MB) were committed to this branch and removed two commits later; a merge commit would land them on main permanently, a squash drops them. .gitignore now covers *.npy/*.npz so it cannot recur.

Scope

Three things, in one branch because they share the README:

  1. Two Binder-runnable reader notebooks + notebooks/viewers.py.
  2. .binder wiring for them.
  3. Archives the nine stale example notebooks off main (~3,832 lines) — custom_aggregations, rasterized_zarr, jupyterhub_example, cryocloud_example, cost_reporting, aoi_mask, shardmap_viewer, tdigest_reader_example, sentinel2_fusion. They were written against older APIs and had drifted; stale documentation is worse than none. Preserved at commit c56221b4 (also on the claude/archive-notebooks-2026-08-26 branch), and the four doc links that pointed at them now use the commit permalink.

Two reader-only demo notebooks, Binder-runnable, plus the .binder and README wiring they need.

notebook what it shows
notebooks/hhdc_viewer.ipynb polygon → MOC → shard → rotatable paired 3-D view (ATL03 + GEDI) → numpy tensors
notebooks/waveform_viewer.ipynb the cell-level join: one GEDI o18 footprint against the 2×2 ATL03 o19 cells beneath it, both from stored digests

Both read the anonymous public source.coop demo stores and import only mortie + moczarr[zagg] — no zagg public API calls, no hand-rolled bit arithmetic. (They are not zagg-free: moczarr imports the t-digest algebra from zagg rather than vendoring it, which is exactly what the [zagg] extra carries. Earlier wording claimed "nothing from zagg itself" and was wrong.) They are split in two because one needs %matplotlib widget for the rotatable 3-D view and the other %matplotlib inline; the two backends collide in a single kernel.

Phases

  • 1 — the two notebooks + viewers.py. Reader-only demos sharing the drawing code, so the cells stay about the read path.
  • 2 — .binder + README wiring. postBuild gets the frontend-extension package a kernel cannot install for itself; README gains the section and both Binder badges.
  • 3 — binned tensor mode in the 3-D view, checked by default, with the exact centroids swept lazily behind it.
  • 4 — archive the nine stale notebooks off main and repoint the docs/ references. This phase is a scope change not described when the PR opened — see the questions below.
  • 5 — fold the adversarial self-review. All fifteen inline threads answered; fourteen resolved, one left standing (the scope point above). Details below.

Dependencies: %pip, not pyproject.toml (espg, 2026-08-26)

Each notebook carries its own %pip install line rather than a new zagg extra, so pyproject.toml is untouched and there is no zagg[demo] -> moczarr[zagg] -> zagg cycle. They therefore run outside Binder unchanged.

One exception lives in .binder/postBuild: ipympl ships a Jupyter frontend extension, which a running kernel cannot pick up from a %pip in a cell. It was two before the review — ipywidgets already arrives with the viz extra, since ipyleaflet>=0.19 requires ipywidgets<9,>=7.6.0, and re-installing it unpinned would one day climb past that cap with a green build.

The moczarr>=0.7 pin is load-bearing. hhdc.block_rank and cell_index land in moczarr 0.7.0 (espg/moczarr#52, espg/moczarr#53). Unpinned, pip resolves 0.6.0 and the notebook dies deep in a cell with an ImportError; pinned, it fails at the install line with something a reader can act on.

The bug this fixes on the way

hhdc_viewer's ancestor hand-rolled the morton bit decode and skipped point_to_area29, so a located companion's level-28/29 digits decoded out of range and rank_to_xy raised rank must lie in [0, 4). _grid_xy is now three library calls — block_rank (which normalizes point words internally), rank_to_rowcol, and the scaling to meters — with no bit arithmetic left in the notebook.

What the review round changed

Nineteen commits, one per finding. The substantive ones:

Correctness

  • grid_xy returned (row, col) while calling it (x, y), so view3d's "east"/"north" labels sat on the wrong axes. It now returns (col, row), agreeing with rank_to_rowcol's documented (row, col) = (y, x), and _binned_pts makes the same swap so both paths still agree. This transposes the rendered scatter relative to every previous run of this branch — intended, but visible, and worth an eyeball on the first real Binder pass.
  • _binned_pts placed points at cell corners while the exact path used cell centres, so toggling the binned box shifted the whole cloud by half a cell — 6.2 m for ATL03's o19 cells, 12.4 m for GEDI's o18, i.e. a different shift per sensor. Unticking that box is the notebook's advertised comparison, so the one thing a reader is invited to do was the thing that exposed an offset that is not real.
  • min(photons, pe) did not rank by "the weaker member": GEDI photoelectrons run two to three orders of magnitude above ATL03 photon counts, so np.minimum collapsed to the ATL03 count on essentially every joint cell. Each side is now scaled by its own maximum over the joint cells before the min, and the prose says what the code does — including that what makes a pick coincident is the joint = (A2 > 0) & (G2 > 0) mask, not the min.
  • _atl03_digests used one parameter, gside, for two different quantities that are equal only for this store pair: GEDI's cells per block edge (2**(18-12)) and ATL03's cells per chunk edge (2**(19-13)). Split, with the chunk side derived from the cell order already on the field's path and block_order + 1 replaced by an explicit atl03_chunk_order.
  • The empty-digest guard in _atl03_digests was defeated three lines later by min() on a zero-size array. An empty cell now draws the GEDI side alone and titles itself accordingly instead of raising inside an Output widget.

Binder survivability (mybinder.org caps the whole container at 2 GB)

  • view3d prefetched 1.28 GiB of tensors — every block, both sensors — to label a dropdown, though the view only ever draws one block. It now reduces each block as it arrives to the drawable cloud (already capped at cap) plus the occupancy count, which is well under 100 MiB and no longer scales with n_bins.
  • hhdc_viewer cell 10 held ~512 MiB of cubes before np.savez_compressed. The cubes are now streamed into an open ZipFile one at a time; an .npz is a zip of .npy members, so np.load reads the result back identically (verified).
  • Both interactive views leaked a matplotlib figure per widget interaction — under %matplotlib widget those are live canvases holding a websocket comm and a mouse handler for the kernel's lifetime. Each view now closes the previous figure (not the current one, which the widget backend still needs to render).

Honesty and prose

  • export asked for resolution=0.5 and got 1 m, and asked for the default 1.0 and got 4 m, because fit="degrade_resolution" is export's default — while printing a sentence that read as though the request had been honoured, under a comment saying # default 64 x 1 m bins. It now compares dz to the requested resolution and says which happened; the fit override is documented.
  • The header claimed "two pip installs" where there is one %pip install line naming five packages.
  • waveform_viewer pointed the reader at 07_minimal, which is untracked staging and not in the tree. Both notebooks now cross-link each other.
  • README said the notebooks import "nothing from zagg itself" (untrue), and repeated the whole reader-only paragraph twice.

Hygiene

  • .gitignore gained notebooks/*.np[yz]. hhdc_viewer cells 9–10 write exactly those files into notebooks/ every time anyone runs it to the end, and nothing stopped a recurrence.
  • The archived-notebook links pointed into claude/archive-notebooks-2026-08-26 — the one namespace reserved for agent working branches and routinely cleaned up. All four now use the commit permalink c56221b4, this PR's own base and the last commit on main that carried them.
  • Both notebooks were red under ruff check (5 errors: two I001, unused matplotlib.pyplot in both, unused numpy in waveform_viewer) and under ruff format, neither of which is excluded for .ipynb in .pre-commit-config.yaml. Both are clean now.

Merge note: the branch history carries two .npy blobs

347e72e8 committed notebooks/atl03_4331422233111.npy (4.0 MB) and notebooks/gedi_4331422233111.npy (2.0 MB); bb797b2a removed them. They are gone from the tree but permanently in this branch's history, so a merge commit would carry them onto main. A squash merge is the remedy — force-pushing the branch to strip them is out of the question. Flagging it here rather than leaving it to be discovered afterwards; the merge strategy is the maintainer's call, not something this PR should decide.

How it was tested

  • _grid_xy pinned against a mz.morton_decimal digit oracle on live ATL03 located words: 0 mismatches, and byte-identical (0.00e+00 m) to the pre-block_rank form.
  • waveform_viewer's full join driven headlessly against the live public store: 64 blocks per sensor, 63 joint, top cell pairing 199 GEDI centroids (12,303 pe) against 97 ATL03 photons, CDFs monotonic and normalized 0→1.
  • The adversarial review executed every code cell of both notebooks in order against the live store (%pip/%matplotlib stripped, Agg backend); the memory and timing figures quoted above are that run's measurements.
  • The streamed .npz was checked to round-trip through np.load with keys, values, and dtype preserved.
  • Local green bar after the fold: pytest -v4891 passed, 44 skipped. ruff check/ruff format clean on all three files this PR touches under notebooks/. pre-commit run --all-files shows no new failures — the remaining ones are pre-existing on main in files this PR does not touch (N818 in src/zagg/registry.py, ruff format on tests/data/benchmark/README.md, ~177 mypy errors under files: src|tests, codespell in src/+docs/specification.md, check-yaml on the CloudFormation !Equals tags).
  • Both notebooks parse, carry cleared outputs, and hold every import in the first code cell.

Questions for review

  1. Still never run top-to-bottom in a fresh kernel on Binder. Every piece is verified headlessly, but the widget backend end-to-end is exactly the claim Binder makes. Two review fixes make this more than a formality: the axis swap transposes the rendered scatter, and the figure-registry bound changes the widget lifecycle. Both are argued from the backend's documented behaviour rather than observed in a live kernel.
  2. Phase 4 is a scope change the PR did not open with. bb797b2a deletes nine notebooks (3,832 lines) from main and rewrites three docs/ pages, under a PR described as "two reader-only demo notebooks plus the wiring they need". The archive itself may well be right — those notebooks had drifted against older APIs — but it is a separate decision from shipping these two, and it is not a review round's to absorb. Either is fine by me: (a) keep it and restate the PR's scope, or (b) split it into its own PR and drop it from here. This is the one review thread left unresolved.
  3. Two postBuild extras may now be orphaned. With jupyterhub_example, shardmap_viewer and cryocloud_example archived, the catalog and viz extras installed by .binder/postBuild have no consumer left in the tree, and their comment still names those notebooks. Removing an extra is a dependency change, so nothing was touched — but note viz is currently what supplies ipywidgets (via ipyleaflet), so dropping it would need the ipympl line to name ipywidgets again, pinned inside <9.
  4. moczarr[zagg] can shadow the Binder checkout. The extra resolves zagg>=0.40, while postBuild's tag-less fallback installs the checkout as 0.0.0+binder, which does not satisfy it — so on a build where no tag is reachable the notebooks' %pip install quietly pulls a released zagg from PyPI over the checked-out one, in shared site-packages, for every notebook in the session. Documented in postBuild rather than fixed, because the obvious fix (drop [zagg] from the notebook install lines) would break the notebooks running outside Binder unchanged.
  5. Deriving the chunk geometry from the store (stored_chunk_spans / iter_populated_chunks) is the better version of the gside split, and would remove the last assumed constant. Not done: moczarr is not installable in this environment, and a wrong guess at a store-addressing API surfaces as silently mis-addressed cells rather than an error. atl03_chunk_order is now the single place to change it.

moczarr 0.7.0 is not on PyPI at the time of writingresolved: 0.7.0 is published (latest on PyPI as of 2026-08-27), so the >=0.7 pin now resolves and the install line no longer fails by design.

Comment thread README.md Outdated
Comment thread .binder/postBuild Outdated
Comment thread .binder/postBuild Outdated
Comment thread notebooks/viewers.py
Comment thread notebooks/viewers.py Outdated
@espg

espg commented Aug 26, 2026

Copy link
Copy Markdown
Member Author

🤖 from Claude (review)

Both notebooks fail ruff — so pre-commit run --all-files (CLAUDE.md §4's local green bar) is not green on this branch.

ruff-pre-commit v0.14.10's ruff-check hook covers jupyter file types, and it is configured here with no .ipynb exclusion (only ^archive/), so the notebooks are linted. Run with the same flags as the PR lint bot (.github/workflows/lint.yml):

$ ruff check --select=E,F,W,I --ignore=E501 notebooks/hhdc_viewer.ipynb notebooks/waveform_viewer.ipynb
I001 Import block is un-sorted or un-formatted   hhdc_viewer.ipynb:cell 2:4:1
F401 `matplotlib.pyplot` imported but unused     hhdc_viewer.ipynb:cell 2:7:29
I001 Import block is un-sorted or un-formatted   waveform_viewer.ipynb:cell 2:4:1
F401 `matplotlib.pyplot` imported but unused     waveform_viewer.ipynb:cell 2:6:29
F401 `numpy` imported but unused                 waveform_viewer.ipynb:cell 2:7:17
Found 5 errors.

All five are real, and the F401s are a direct consequence of the viewers.py refactor: the drawing moved out, so plt is now dead in both notebooks and np is dead in waveform_viewer (nothing in its cells 3/5/7 touches numpy — cell 5 uses only ndarray methods). Note hhdc_viewer still uses np for np.save/np.savez_compressed, so only plt goes there.

The I001s are moczarr/mortie being sorted as third-party while viewers sorts as first-party — ruff's autofix produces the right grouping.

(notebooks/viewers.py itself is clean: ruff check --select=E,F,W,I,N and ruff format --check both pass.)

@espg

espg commented Aug 26, 2026

Copy link
Copy Markdown
Member Author

🤖 from Claude (review)

notebooks/hhdc_viewer.ipynb — three findings from a headless top-to-bottom run.

Executed every code cell in order against the live public store (%pip/%matplotlib lines stripped, Agg backend, viewers.py at 1b96153). All cells completed; the numbers below are that run's actual output.


(1) Cell 9 silently reports a z resolution it did not get, and its inline comment is wrong. Measured output:

gedi: (64, 64, 128) tensor, z = -38.0 m + bin * 1 m
atl03: (128, 128, 64) tensor, z = -68.0 m + bin * 4 m

The GEDI call asks for resolution=0.5 and gets 1 m; the ATL03 call asks for the default 1.0 and gets 4 m. fit="degrade_resolution" is the export default, so this is working as designed — but the print interpolates dz into a sentence that reads like the request was honoured, and the trailing comment on the ATL03 call says something that is simply false for this store:

atl03, meta = export(view, "atl03")           # default 64 x 1 m bins

Cell 10 gets this exactly right two cells later ("48 at the requested 1 m bins; 16 degraded to fit their relief in 128 bins"), which makes the gap more conspicuous, not less. export's docstring likewise says "n_bins/resolution set the vertical shape" without mentioning that fit can override resolution. Suggest having export compare dz to resolution and say so, and fixing the comment.

(2) Cell 10 holds ~512 MiB of cubes in memory at once, before compressing. Measured: 64 blocks × (128, 128, 128) uint32 = 8.0 MiB each, all retained in cubes and then handed to np.savez_compressed in one call. mybinder.org's default cap is 2 GB for the whole container, and the kernel is already holding two open leaf stores plus view3d's prefetched centroids (1.27 M ATL03 + 1.76 M GEDI). It survived here, but it is the one cell in the PR that could OOM a Binder session, and nothing warns the reader. Streaming block-by-block into an open ZipFile, or capping the number of cubes kept, would remove the risk. (The rest of the science in that cell checks out: an o12 block really is 2**(19-12) = 128 ATL03 cells a side, and the shapes came back (128, 128, 128).)

(3) The shard is swept twice per sensor over the network. Cell 5's open_shard runs a full read_ragged sweep purely to print a centroid count and throws the data away; cell 7's view3dviewers.load then re-reads the same two arrays end to end. Measured cost of the discarded half:

atl03  open+sweep  14.7s (+  48 MB) — 1,266,765 centroids
gedi   open+sweep  13.4s (+   1 MB) — 1,764,480 centroids

so ~28 s of the Binder demo is spent reading bytes that are immediately dropped and re-fetched. Either have open_shard return what it read and hand it to view3d, or drop the sweep and report what open_ragged already knows.


Prose: the header says "in five functions, two pip installs, and zero credentials". There is one %pip install line, and it names five packages. The closing cell's "Five functions, two libraries" is the claim that actually holds — the header looks like it drifted from an ancestor. ("Reader-side — mortie for the geometry, moczarr for the store" is also not quite right; see the README thread.)

@espg

espg commented Aug 26, 2026

Copy link
Copy Markdown
Member Author

🤖 from Claude (review)

notebooks/waveform_viewer.ipynb — two findings.

Ran every code cell in order against the live public store (%pip/%matplotlib stripped, Agg backend). It works, and reproduces the PR body's numbers: 4 candidate shards, 64 blocks per sensor, 63 joint, dropdown built and populated.


(1) Markdown cell 0 points the reader at a notebook that does not exist in this repo.

It is a separate notebook from 07_minimal on purpose: that one needs %matplotlib widget

07_minimal is untracked staging (demo/07_minimal.ipynb, demo/07b_minimal.ipynb) — it is not in the tree, and this is the only reference to it anywhere:

$ git grep -n "07_minimal"
notebooks/waveform_viewer.ipynb:15:    "It is a separate notebook from `07_minimal` on purpose: that one needs\n",

A Binder reader clicking through has no way to find it. The sibling shipped by this PR is hhdc_viewer.ipynb, and naming it also makes the pair discoverable in both directions (hhdc_viewer's header never mentions waveform_viewer either).

(2) min(photons, pe) does not rank by "the weaker member" — the two units are ~500× apart. Cell 6's prose says:

The slider ranks joint cells by the weaker member — min(photons, pe) — so early picks are genuinely coincident rather than one-sided.

but viewers.py:417 takes np.minimum across incommensurate units. From cell 5's own measured output on this shard:

  4331422233341    90 joint o18 cells       1,871 atl03 photons     886,872 gedi pe
  4331422233344    84 joint o18 cells       2,098 atl03 photons     582,428 gedi pe
  4331422233313    71 joint o18 cells       2,884 atl03 photons     402,081 gedi pe

GEDI photoelectrons run two to three orders of magnitude above ATL03 photon counts, so np.minimum(A2, G2) reduces to A2 on essentially every joint cell and the slider is really "rank by ATL03 photon count". What actually delivers the "genuinely coincident rather than one-sided" property is the joint = (A2 > 0) & (G2 > 0) mask one line earlier — which is a good property, just not the one the sentence attributes to the min.

Either normalize each side before taking the min (e.g. rank on min(A2/A2.max(), G2/G2.max()) over the joint cells) so it means what it says, or keep the current ordering and describe it accurately: "cells both sensors populate, ordered by ATL03 photon count".

Comment thread notebooks/viewers.py
Comment thread README.md Outdated
Comment thread notebooks/viewers.py Outdated
Comment thread README.md Outdated
Comment thread notebooks/viewers.py Outdated
Comment thread README.md Outdated

They share `notebooks/viewers.py`, which holds the drawing so the notebooks stay about the read path. The two are split because one needs `%matplotlib widget` for its rotatable 3-D view and the other `%matplotlib inline`; the backends collide in a single kernel.

**Archived notebooks.** The earlier examples (`custom_aggregations`, `rasterized_zarr`, `jupyterhub_example`, `cryocloud_example`, `cost_reporting`, `aoi_mask`, `shardmap_viewer`, `tdigest_reader_example`, `sentinel2_fusion`) were written against older APIs and had drifted out of date — stale documentation being worse than none. They are preserved on the [`claude/archive-notebooks-2026-08-26`](https://github.com/englacial/zagg/tree/claude/archive-notebooks-2026-08-26/notebooks) branch and removed from `main`.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤖 from Claude (review)

Four permanent doc links now point into a claude/* working branch, which is the one namespace guaranteed to be unstable.

This line plus docs/aoi_mask.md:201, docs/quickstart.md:67, and docs/quickstart.md:205,215 now resolve through claude/archive-notebooks-2026-08-26. The branch does exist right now —

$ git ls-remote --heads origin claude/archive-notebooks-2026-08-26
c56221b4396ce386836d313137c1f383eb1240d8  refs/heads/claude/archive-notebooks-2026-08-26

— but claude/* is exactly the namespace CLAUDE.md §2 reserves for agent working branches and that routine runs create, push to, and (post-merge) clean up. Published documentation on main pointing at it will break silently the first time someone tidies stale agent branches, and nothing in the branch name says "do not delete me".

Two stable alternatives, either fine:

  • A commit permalink. c56221b4 is the last commit where those notebooks lived on main, and GitHub keeps blob URLs at an explicit sha forever, independent of any branch: https://github.com/englacial/zagg/blob/c56221b4/notebooks/aoi_mask.ipynb. Zero new refs to maintain.
  • A tag (e.g. notebooks-archive-2026-08-26) — but note tag creation in this repo is a production trigger (publish.yml fires on *.*.*), so it must not look like a version, and §1 puts tag pushes out of an agent's hands anyway.

Separately: this is a scope change the PR body does not describe. The body opens "Two reader-only demo notebooks … plus the .binder and README wiring they need"; bb797b2 additionally deletes nine notebooks (3,832 lines) from main and rewrites three docs/ pages. Under §2 that belongs in the description with its own phase entry — and under §6 a scope change of this size is an @espg decision rather than something a review round should quietly absorb.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤖 from Claude

Split this one: the link half is folded, the scope half is not mine to take and stays open.

Folded (62ca510) — your first alternative, the commit permalink. All four links now resolve through c56221b4 instead of claude/archive-notebooks-2026-08-26:

  • README.md:147, docs/aoi_mask.md:201, docs/quickstart.md:67, docs/quickstart.md:205 and :215.

Verified c56221b4 actually carries them before repointing — git ls-tree c56221b4 notebooks/ lists all nine (aoi_mask, cost_reporting, cryocloud_example, custom_aggregations, jupyterhub_example, rasterized_zarr, sentinel2_fusion, shardmap_viewer, tdigest_reader_example). It is this PR's own base, i.e. the last commit on main that had them. The README prose changed too, since it said "branch" while pointing at a sha:

They are removed from main and preserved at c56221b4, the last commit on main that carried them — a commit permalink rather than a branch link, because GitHub keeps a blob URL at an explicit sha forever and there is no ref anyone has to remember not to delete.

I did not create a tag; agents do not push tags here, and a *.*.* tag triggers publish.yml.

Left standing — the scope point. You are right that bb797b2a deletes nine notebooks and rewrites three docs/ pages while the PR body still opens "two reader-only demo notebooks … plus the .binder and README wiring they need". That is a scope change, and a review round is not what should absorb it. I have described it in the PR body under the phases checklist and raised it as an explicit question rather than acting further on it — reverting the archive, or keeping it and restating the PR's scope, is a call for the maintainer, not for this fold.

Leaving this thread unresolved for that reason.


Generated by Claude Code

@espg

espg commented Aug 26, 2026

Copy link
Copy Markdown
Member Author

🤖 from Claude (review)

Two .npy blobs (6 MB) are in this branch's history, and .gitignore does not stop it happening again.

347e72e committed them and bb797b2 removed them, so both are permanently in the PR's commit history:

$ git log --oneline --diff-filter=A -- notebooks/atl03_4331422233111.npy notebooks/gedi_4331422233111.npy
347e72e8 viewers: binned tensor mode in the 3-D view, checked by default

$ git show bb797b2a --stat | grep npy
 notebooks/atl03_4331422233111.npy      | Bin 4194432 -> 0 bytes
 notebooks/gedi_4331422233111.npy       | Bin 2097280 -> 0 bytes

Catching them was right — CLAUDE.md §3 is explicit about generated artifacts and large binaries. Two things still open:

(1) They will land on main unless this merges as a squash. A merge commit brings the whole branch history, blobs included; §1 rules out force-pushing the branch to strip them, so the squash is the only remedy and it is an @espg call at merge time. Worth a line in the PR body so it is not discovered afterwards.

(2) Nothing prevents a recurrence. .gitignore has no *.npy / *.npz rule:

$ grep -nE 'npy|npz' .gitignore
(no matches)

and hhdc_viewer.ipynb cells 9–10 write exactly these files into notebooks/, by construction, every time anyone runs the notebook to the end:

np.save(f"gedi_{meta['block']}.npy", gedi)
np.save(f"atl03_{meta['block']}.npy", atl03)
np.savez_compressed(f"atl03_cubes_{shards[0]}.npz", ...)

The .npz is the bigger one — 64 compressed (128,128,128) uint32 cubes. A notebooks/*.np[yz] entry (or a global *.npy / *.npz) closes it. Writing the exports to a gitignored subdirectory, or to tempfile.gettempdir(), would be the alternative — but the notebook's point is partly that the tensors land somewhere a reader can pick them up, so the ignore rule is probably the better fit.

Comment thread notebooks/viewers.py Outdated
Comment thread notebooks/viewers.py Outdated
Comment thread notebooks/viewers.py Outdated
Comment thread notebooks/viewers.py
@espg

espg commented Aug 26, 2026

Copy link
Copy Markdown
Member Author

🤖 from Claude (review)

Review bookkeeping — the head moved mid-review, and one earlier point is now superseded.

This round started on 1b961539 and the branch advanced to bb797b2a (347e72e8 binned-tensor mode, bb797b2a archive + .npy removal) while it was running. Re-verified everything against bb797b2a:

Still valid, unchanged code:

  • ruff / pre-commit — same 5 errors at bb797b2a (hhdc_viewer cell 2 plt; waveform_viewer cell 2 plt + np; two I001).
  • hhdc_viewer.ipynb and waveform_viewer.ipynb are byte-identical to 1b961539, so every notebook finding stands: cell 9's silent 4 m degradation, cell 10's 512 MiB, the "two pip installs" header, the dangling 07_minimal reference, the min(photons, pe) ranking.
  • .binder/postBuild is untouched by both new commits, so both threads there anchor correctly.
  • viewers.py: grid_xy (orientation), _atl03_digests (gside), waveform_view (figure accumulation, binw refetch, empty-digest min) are all unchanged.

Superseded — please disregard point (3) of the hhdc_viewer comment. I wrote that cell 5's discarded sweep is immediately repeated by view3dload. 347e72e8 removed that eager load and made it lazy behind _exact(), so the double ragged sweep only happens now if the reader unticks "binned". The half that still stands is that cell 5's 28 s sweep is thrown away outright — and 347e72e8 replaced it with something more expensive, which is its own thread: 88 s and 1.28 GiB of tensors read up front to label a dropdown.

Three threads were mis-anchored — posted against 1b961539 using line numbers already read from bb797b2a. Each carries a reply pointing at its correctly-anchored replacement; act on the replacements:

stale replacement subject
r3866658469 r3866764104 grid_xy returns row/col, not x/y
r3866660963 r3866764267 gside conflates two different 64s
r3866673708 r3866764435 figures never closed

Two new findings came in with 347e72e8 / bb797b2a and are on their own threads: the 1.28 GiB tensor prefetch (the one I would fix first), and the half-cell xy offset between binned and exact mode.

Not findings, checked and clean: SIDE is right (sqrt(4π/(12·4¹²))·R = 1591.8 m); an o12 block really is 128 ATL03 o19 cells a side and the cubes came back (128, 128, 128); A2.sum(axis=2).reshape(64,2,64,2).sum(axis=(1,3)) is a correct o19→o18 fold given the nested-rank row/col mapping; cell_index(chunk, …) indexing into generate_morton_children is correct (mortie generates children in nested order); dividing cdf_from_tdigest by the total weight is the right normalization (it returns cumulative weight with endpoints 0 and total); both Binder badge URLs are well-formed and point at paths this PR creates; notebooks/viewers.py passes ruff check --select=E,F,W,I,N and ruff format --check; and viewers.py carries no leftover notebook-only globals — display is imported from IPython.display, interactive_output calls its callback once on creation so view.block is set before cell 9 runs, and the cdf_from_tdigest import is genuinely function-local. waveform_view's draw callback builds a real 3-axis figure (verified headlessly at pair=0/nth=0/binw=1.0, pair=0/nth=5/binw=0.0, and pair=62/nth=39), so nothing silently no-ops.

@espg espg added the implement label Aug 27, 2026
@espg

espg commented Aug 28, 2026

Copy link
Copy Markdown
Member Author

🤖 from Claude (review)

Verdict: the ATL03 × GEDI waveform join is spatially sound. Scoped review of notebooks/viewers.py (waveform_view / _atl03_digests / _binned_pts) and notebooks/waveform_viewer.ipynb's pairing cell, hunting for another instance of the b9d7631e failure class (silent reorder / mis-index that a consumer assumes is aligned). Every check below was run against the live anonymous store (shards[0] = 4331422233, block 4331422233341, 90 joint cells, moczarr 0.7.0 / mortie 0.9.11 / zagg 0.52.0), not read off the page.

(1) The o19 → o18 fold groups the right quads. at.sum(axis=2).reshape(GSIDE, 2, GSIDE, 2).sum(axis=(1, 3)) on a 128×128 tensor whose values encode row*1000 + col reproduces sum over (2r+dr, 2c+dc) exactly at every probed (r, c). The plausible-looking alternative reshape(2, GSIDE, 2, GSIDE).sum(axis=(0, 2)) does not equal it, so the written order is load-bearing and correct.

(2) The fold's (r, c) is the same ground as GEDI's (r, c). Independent ground-truth maps were built for the top block straight from read_ragged — cell word → block_rankrank_to_rowcol → accumulate stored weight — and compared against the read_tensors tensors under identity and seven index perturbations:

variant ATL03 o19 128×128 (corr / occupancy Jaccard) GEDI o18 64×64
identity +0.9911 / 1.0000 +1.0000 / 1.0000
transpose −0.0230 / 0.0423 −0.0013 / 0.0658
flipud −0.0111 / 0.0487 −0.0143 / 0.0505
fliplr −0.0078 / 0.0451 +0.0210 / 0.0717
roll row ±1 +0.16 / 0.299 −0.031 / 0.061
roll col ±1 +0.31 / 0.382 −0.061 / 0.030

Occupancy matches cell-for-cell (1529/1529 ATL03, 478/478 GEDI). Both tensors are indexed by the nested-rank deinterleave at their own depth, and rank_to_rowcol satisfies parent(row, col) == (row//2, col//2) over all 4**7 o19 ranks — which is exactly what makes A2[r, c] and G2[r, c] the same patch of ground, same origin corner, no transpose.

(3) _atl03_digests(w, r, c) reads exactly the four cells the fold summed. For the widget's 10 top-ranked joint cells on the top block, the total weight of the digests actually fetched through generate_morton_childrenrowcol_to_rankcell_indexread_cell was compared against the ground-truth map:

nth   r   c   read_wt  AW[2r,2c]  AW[2c,2r] AW[2r+1..]       AW.T  A2(tensor)
  0  46  59      84.0       84.0        0.0       32.0        0.0        85.0
  1  29   5      57.0       57.0        0.0      100.0        0.0        56.0
  3  63  49      52.0       52.0       29.0       41.0       29.0        48.0
  4  27   2      97.0       97.0        0.0      108.0        0.0        92.0
  ...
identity 10/10 exact; transposed-rc 0/10; off-by-one row 0/10; AW.T 0/10

Identity is the outright winner — exact, to the photon, on every pick — and every transposed/offset variant fails. The residual A2 / read_wt ≈ 0.931 is the tensor's finite z window, as the docstring says (see note (4) in the separate comment).

GEDI's own read path is exact too: read_cell(cell_index(store, field, block_word, r, c)) matched GW[r, c] to the milli-pe on 5 probes (19728.420 / 14936.512 / 13923.032 / 12434.455 / 12302.557), and GW[c, r] was 0.000 on all five.

(4) Chunk addressing is right. On the top block, kids[rowcol_to_rank(kr, kc, depth=1)] decodes back to (kr, kc) for all four quadrants:

(row=0,col=0) -> 43314222333411 -> rowcol 0,0
(row=0,col=1) -> 43314222333412 -> rowcol 0,1
(row=1,col=0) -> 43314222333413 -> rowcol 1,0
(row=1,col=1) -> 43314222333414 -> rowcol 1,1

No off-by-one in morton child order. (Separate comment on what the bare except does and does not protect.)

(5) Sortedness holds on every path into cdf_from_tdigest. b9d7631e's sort fixes the ATL03 merge; the GEDI side is a single stored digest. Swept the whole leaf: 0 of 64,499 ATL03 and 0 of 18,774 GEDI stored digests have a decreasing mean anywhere (369 and 1 respectively have tied means, which np.interp handles as a vertical step — correct behaviour, not a defect). cdf_from_tdigest is the only interp-like consumer in the file — _mixture is a weighted Gaussian sum and np.histogram bins by value, neither of which cares about order.

(6) The rank normalization does what it says. Both maxima are taken over [joint], not over all cells. Division cannot blow up: max(int(...max()), 1) floors the divisor at 1, and pairs only holds blocks where joint.any(), so a2[joint].max() is never called on an empty selection. Both tensors are uint32 (folded A2 uint64), so the int() cast is lossless here (measured truncation 0.0 on both). Over 200 randomized joint layouts plus the live block, np.argsort(rank.ravel())[::-1] with min(nth, joint.sum() - 1) always lands on a joint cell — non-joint cells are forced to rank 0 by * joint, joint cells are strictly positive because A2 > 0 and G2 > 0, so the joint cells occupy exactly the first joint.sum() slots of the descending order and the clamp cannot walk past them. np.unravel_index recovers the cell the score belongs to (rank[r, c] == rank.ravel()[idx] on every trial).

Reproduction scripts are throwaway; the numbers above are all from live reads of s3://us-west-2.opendata.source.coop/englacial/zagg/demo/{atl03_tdigest_o9,gedi_flux_o9}.zarr.

@espg

espg commented Aug 28, 2026

Copy link
Copy Markdown
Member Author

🤖 from Claude (review)

Finding (low, hardening): waveform_view's gside parameter is documented, passed by the notebook, and never used — while the 2× coupling it describes is hardcoded in three places with nothing asserting it.

notebooks/viewers.py, waveform_view:

def waveform_view(
    stores, fields, blocks, pairs, shard, gside: int = 64, atl03_chunk_order: int = 13
):
    """...
    ``gside`` is GEDI's cells across an o12 BLOCK edge. ...
    """

AST walk of the function body (ast.Name ids inside the waveform_view FunctionDef):

gside referenced in waveform_view body: False
atl03_chunk_order referenced: True

notebooks/waveform_viewer.ipynb passes it explicitly (waveform_view(stores, FIELDS, blocks, pairs, shards[0], gside=GSIDE)), and _atl03_digests's docstring builds a careful distinction around it — "That is deliberately NOT the caller's gside... Sharing one name would make a change to either sensor's cell order... address the wrong cell silently rather than raise" — so a reader reasonably concludes the o18 grid width is being honoured somewhere. It is not: r/c come from rank.shape, which comes from the caller's G2.

Benign today, but it is the one piece of geometry in this join that is asserted nowhere. The ATL03:GEDI 2× ratio is hardcoded three times:

  • the notebook: .reshape(GSIDE, 2, GSIDE, 2).sum(axis=(1, 3))
  • _atl03_digests: rr, cc = 2 * r + dr, 2 * c + dc over dr, dc in (0, 1)
  • the title string: f"vs ATL03 {len(adigest)} (2×2 @o19)"

waveform_view already derives acell_order = int(fields["atl03"].split("/", 1)[0]) from the field path. A store pair at o19/o17 (or o18/o17) would make the notebook's reshape and _atl03_digests's 2×2 walk disagree about which ground each (r, c) covers, and — as the companion comment on the except shows — that disagreement produces plausible wrong numbers, not an exception.

Two concrete options:

  1. Use it. Derive gcell_order the same way acell_order is derived (int(fields["gedi"].split("/", 1)[0])), assert 2 ** (acell_order - gcell_order) == 2 (or generalize the walk to range(step)), and drop the unused gside parameter — the notebook's GSIDE then only drives the reshape, where it is genuinely used.
  2. Drop it. Delete the parameter and the notebook's gside=GSIDE, and reword the _atl03_digests docstring paragraph so it no longer contrasts chunk_side against a parameter that does not exist.

Related nit in the same spirit: ASIDE in waveform_viewer.ipynb cell 1 (ASIDE, GSIDE = 128, 64) is defined and never referenced (1 occurrence in the whole notebook's code — its own definition), while GSIDE is used 4 times.

Not a correctness bug — flagged because it is the only un-asserted assumption left in an otherwise verified join.

@espg

espg commented Aug 28, 2026

Copy link
Copy Markdown
Member Author

🤖 from Claude (review)

Assessment (low, requested explicitly): the bare except (KeyError, ValueError): pass in _atl03_digests does not hide addressing errors — but only because an addressing error would not raise at all.

notebooks/viewers.py, _atl03_digests:

except (KeyError, ValueError):
    pass  # that quarter holds no photons

I probed what a wrong quadrant actually does. Taking an occupied o19 cell that lives in quadrant (1, 1) of block 4331422233341 and reading it through the correct child versus child 0:

cell (64,80) via correct quadrant 43314222333414: 3 centroids, wt=3.0  (truth AW=3.0)
cell (64,80) via wrong   quadrant 43314222333411: 13 centroids, wt=13.0 (truth AW=3.0)

No exception. A wrong-block chunk word behaved the same way — cell_index + read_cell returned an array rather than raising. So the failure mode this guard looks like it is protecting against (an off-by-one in morton child order silently reading a neighbouring quadrant) would sail straight past it as a plausible-looking digest of the wrong ground, exactly the b9d7631e shape.

The guard's stated purpose is sound and it fires legitimately and often: over the 41 picks the nth slider can reach on the top block, 69 of 164 quarters are genuinely empty in the ground-truth map, and 0 picks have all four empty (so the "no ATL03 digest under this cell" banner is a real but unexercised path on this AOI). Addressing itself is verified exact — see the companion verdict comment, 10/10 to the photon.

Two things would make the guard mean what its comment says, cheaply and with no extra I/O:

  1. Assert the decoded chunk word rather than trusting the index arithmetic:

    chunk = int(kids[rowcol_to_rank(rr // chunk_side, cc // chunk_side, depth=depth)])
    krank, korder = block_rank(np.array([chunk], np.uint64), block_order)
    assert rank_to_rowcol(krank, chunk_order - block_order) == (rr // chunk_side, cc // chunk_side)

    (or the cheaper equivalent, since generate_morton_children is verified nested-rank ordered: nothing more than a one-line comment recording that kids is indexed by nested rank, which is the property the whole quadrant map rests on.)

  2. Narrow the comment to what it can prove: the except catches misses, and cannot distinguish "no photons here" from "asked the wrong place" — which is worth saying out loud in a demo notebook whose whole point is that the join is trustworthy.

No change is strictly required; posting because the question was asked directly and the answer ("the except is not the risk; the silent non-raising read is") is worth having on the record.

@espg

espg commented Aug 28, 2026

Copy link
Copy Markdown
Member Author

🤖 from Claude (review)

Finding (low, labelling): the block table's "atl03 photons" is z-window-clipped and reads ~9% low against the raw digests it sits beside.

notebooks/waveform_viewer.ipynb, pairing cell:

print(
    f"  {mz.morton_decimal(w)}  {int(j.sum()):4,} joint o18 cells   "
    f"{int(A2[j].sum()):9,} atl03 photons   {int(G2[j].sum()):9,} gedi pe"
)

A2 is a sum over the read_tensors tensor, so it counts only what falls inside that block's fitted z window. Measured on block 4331422233341, comparing the tensor against the same block's stored digests read straight through read_ragged:

  • ATL03: tensor.sum() / digest_weight.sum() = 0.9067 for the whole block; 0.9311 mean over the widget's 10 top-ranked joint cells (per-cell A2[r, c] vs the four digests actually read, e.g. 85 vs 84, 56 vs 57, 46 vs 49, 48 vs 52, 92 vs 97, 43 vs 54).
  • GEDI: 0.9999, so this is an ATL03-side effect.

Cause is visible in the fitted windows — the two sensors get different offsets from fit="degrade_resolution" at the same n_bins=256, resolution=1.0: ATL03 (offset, gain) = (-68.0, 1.0)[-68, 188], GEDI (-26.0, 1.0)[-26, 230]. ATL03's noise/blunder photons outside [-68, 188] are dropped from the count. This is by design and the _atl03_digests docstring already anticipates the tensor/digest divergence ("a finite z window, possibly degraded"), but the printed label says "atl03 photons" without qualification, next to a "gedi pe" number that is not clipped, in a cell whose stated purpose is to let a reader judge which blocks have real data on both sides.

It does not affect spatial coherence — the ranking and the cell pick were verified exact — and the ~7-9% shortfall is uniform enough that block ordering is unchanged. Cheapest fix is one word in the f-string: atl03 photons (in-window), or a sentence in the markdown cell above noting that both counts are tensor counts and so bounded by each sensor's own fitted z window. Flagging it because "85 photons" printed against a panel that then plots 84 is the kind of small unexplained gap that costs a reader time.

@espg espg left a comment

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤖 from Claude (review)

Adversarial review of the newest phase only (d217dc8d..0d08fd05). I verified the numerics against the live demo stores with a working moczarr 0.7.0 / mortie 0.9.11 / zagg 0.52.0 env, on the block the notebook actually exports (4331422233444 — confirmed to be view3d's joint[0]).

Verified correct, no comment needed:

  • block_rank(...) >> 2*(29 - order) really is word truncation: ((a29 >> sh) << sh) | order then block_rank gives a bit-identical rank at o19/o22/o24/o27 on all 22,252 located words of the block.
  • Chip naming and decomposition: tiles[rowcol_to_rank(row // side, col // side, depth)] equals the photon's o15 ancestor word for every photon, and (row % side, col % side) equals rank_to_rowcol(block_rank(word, 15) >> 2*(29-22), 7) for every photon.
  • _coverage's reshape(side, k, side, k).any(axis=(1,3)) reproduces the word-derived o18 coverage exactly (702 == 702 parents).
  • registered_pair's placement matches read_tensors independently: atl03 cube.sum(2) > 0 == mask == 2 (1,514 == 1,514), and gedi == np.kron(mask == 2, ones((2,2))) (804 == 804). Cube totals: atl03 22,248 / 22,252 stored; gedi ratio exactly 4.0000005.
  • The mask == 2 fix: astype(bool) gives 1,973 vs 1,514 (+30.3%) and 431 vs 201 (+114.4%) — the commit message's numbers are right to the decimal. No astype(bool) occupancy use remains, and mask == 2 is correct in the 2-state regime too, so no has_exact_occupancy guard is needed.
  • No silent weight loss to np.rint: every ATL03 centroid in this store has weight exactly 1.0 (22,252 centroids for 22,252 photons — nothing merged).
  • chunk_z_range unpack order, and peak RSS ≤ 230 MiB for voxel_chips at o19/o22/o24 and for registered_pair — comfortably inside Binder's 2 GB.

Six findings below, ranked.

Comment thread notebooks/hhdc_viewer.ipynb Outdated
Comment on lines +259 to +268
" normalizes each centroid's order-29 point word and returns its block-local\n",
" nested rank, and shifting that rank right by `2 * (29 - order)` is the same\n",
" word truncated to `order` -- nested ranks are hierarchical, so no re-decode.\n",
" That is the whole reason the cube can be finer than the o19 cells.\n",
"\n",
" The z bin is the cell edge, so a voxel is a cube. Each chip gets its own\n",
" `z0` (the floor of its own minimum) rather than one window for the block:\n",
" at `n_bins` bins a chip spans `side * dz` metres of z, and a block's full\n",
" relief rarely fits in that, so a shared floor would spend most bins empty.\n",
" `z0` is saved per chip, so absolute elevation is always recoverable.\n",

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤖 from Claude (review)

z0 is not actually saved — the .npz has no z datum at all.

The docstring promises "z0 is saved per chip, so absolute elevation is always recoverable", and the comment ten lines down invites the reader to np.load it back ("an .npz is a zip of .npy members, so np.load reads this back unchanged"). But the only things writestr ever puts in the archive are the chip arrays:

zf.writestr(f"{name}.npy", buf.getvalue())
kept[name] = {"z0": float(z0), "voxels": int(chip.astype(bool).sum())}

z0 goes into the in-memory kept dict only. Ran it and checked the artifact:

>>> np.load('atl03_o22_chips_4331422233444.npz').files
['4331422233444111', '4331422233444112', ...]   # 40 members, all chips
>>> [k for k in f.files if not k[0].isdigit()]
[]

So the file on disk carries no z0 and no dz. xy is recoverable (the member names are morton ids), z is not: reopen the file in a fresh kernel and the third axis is an unlabelled bin index. The notebook binds the manifest to manifest in the next cell and nothing persists it.

registered_pair gets this right one cell later — np.savez_compressed(path, **cubes, z0=z0, dz=dz, order=order) — so the fix is the same move: write a z0 member (one float per chip name, or a parallel z0 array plus a scalar dz) into the same zip.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤖 from Claude

Fixed in 7141cf4, and verified by running it rather than by reading it.

voxel_chips now accumulates meta[name] = {"z0": ..., "dz": ..., "written": ...} and writes the manifest into the same zip as a meta.json member:

zf.writestr("meta.json", json.dumps(meta))  # z0 per chip; absolute z is recoverable

Ran the current cell against the live store on the block the notebook now lands on (4331422411132), at both orders, and inspected the artifact:

o22: npz members: 47 | non-chip members: ['meta.json']
     meta.json parsed: 46 entries; sample 4331422411132111 -> {'z0': -20.0, 'dz': 1.5543972149675078, 'written': 3459}
     z0 recoverable from file alone: True
     meta keys == chip members: True
o24: npz members: 481 | non-chip members: ['meta.json']
     meta.json parsed: 480 entries; sample 433142241113211111 -> {'z0': 0.0, 'dz': 0.38859930374187696, 'written': 290}
     z0 recoverable from file alone: True
     meta keys == chip members: True

np.load returns the member as raw bytes (it is not .npy-prefixed, so the loader hands it back verbatim), json.loads parses it, and every chip name in the archive has an entry carrying both z0 and dz — so the third axis is absolute in a fresh kernel with no manifest in scope. The docstring's promise now holds against the file.

Comment thread notebooks/hhdc_viewer.ipynb Outdated
Comment on lines +227 to +231
"| order | voxel | chip covers | xy fill in a 128-chip | full block, dense |\n",
"|---|---|---|---|---|\n",
"| o19 | 12.4 m | 1592 m | 9.2% | 0.01 GiB |\n",
"| **o22** | **1.554 m** | **199 m** | **4.1%** | 0.50 GiB |\n",
"| o24 | 0.389 m | 49.7 m | 0.7% | 8.00 GiB |\n",

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤖 from Claude (review)

The o22 and o24 fill figures don't reproduce; they understate the sparsity by ~3x and ~2x.

Measured with this notebook's own code path on the block it actually exports (4331422233444, which is view3d's joint[0]), using block_rank(locs, 12) >> 2*(29 - order)rank_to_rowcol and counting distinct (row % 128, col % 128) per non-empty 128-chip:

order chips written occupied columns mean fill median max photon-weighted
o19 1 1,514 9.24% 9.24% 9.24% 9.24%
o22 40 9,350 1.43% 0.89% 4.32% 2.58%
o24 298 15,651 0.32% 0.24% 1.43% 0.57%

The o19 row lands on the table's 9.2% exactly, which pins the intended metric as per-chip column fill — and under that same metric o22 is 1.43%, not 4.1%, and o24 is 0.32%, not 0.7%.

I tried to disprove this under every other reading I could construct, and none produces 4.1 / 0.7 either:

  • filled voxels / 128² (from the function's own manifest): 24.05% / 2.35% / 0.34% — and this one breaks the o19 anchor.
  • fill over the whole block rather than per chip: 9.24% / 0.89% / 0.09%.
  • the same statistic over the whole shard (1,266,765 photons, all 64 blocks): mean 6.15% / 0.95% / 0.26%, median 5.72% / 0.79% / 0.19%, max 14.39% / 12.32% / 3.35%.

This matters because the paragraph right under the table is an argument from these numbers ("read the fill before you ask for it … a finer grid buys detail along-track and empty voxels everywhere else") — the real numbers make that argument stronger, and a public demo that overstates o22 fill by 3x is the kind of thing a reader will check.

Everything else in the table is right: 12.435 / 1.554 / 0.389 m voxels, 1592 / 199 / 49.7 m chip extents, 0.01 / 0.50 / 8.00 GiB dense blocks, 8x8 = 64 and 32x32 = 1,024 chips (measured 40 and 298 written).

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤖 from Claude

Fixed in 7141cf4 — but by deletion, not by correcting the figures, so worth saying plainly: the fill column is gone rather than restated at 1.43% / 0.32%.

The condense pass replaced the whole table with two sentences in the export markdown:

Default o22: 1.554 m voxels, z binned to match, the block emitted as 8×8 isotropic 128³ chips (199 m a side), empty chips skipped. order=24 gives 0.389 m voxels in 49.7 m chips — many more, mostly empty.

I checked that no wrong number survives anywhere else. Scanning every markdown cell of the notebook at HEAD for percent-signs and table pipes, the only remaining pipes are in the coverage prose (cell 2) and there is no % claim left in any markdown cell4.1, 0.7% and 9.2% are all absent from the tree. The 4.1% / 0.7% rows last appear at 0d08fd0.

What replaced them makes no unsupported numeric claim: 1.554 m / 199 m / 0.389 m / 49.7 m / 8×8 are the ones you had already confirmed correct, and the sparsity argument is now carried qualitatively ("many more, mostly empty") rather than by a measured fraction. Fewer numbers than the finding asked for, but none of them wrong.

Comment thread notebooks/viewers.py Outdated
Comment on lines +47 to +49

#: What one unit of a sensor's stored WEIGHT is. ATL03 counts photons; GEDI's
#: flux is scaled to photoelectrons. They are not the same quantity and run two

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤖 from Claude (review)

"scaled to photoelectrons" contradicts the store's own attrs and this PR's previous commit.

The demo store's 18/rx_flux declares the gain as an explicit placeholder:

"gain": {"name": "unit", "version": "gedi01b-v002-placeholder", "value": 1.0}

(read live from s3://us-west-2.opendata.source.coop/englacial/zagg/demo/gedi_flux_o9.zarr, leaf 4331422233). So the weights are background-subtracted waveform counts at unit gain — not scaled to photoelectrons.

The sibling notebook says exactly that, and it was written in the commit immediately before this one (d217dc8, "cite the GEDI gain"):

this store ships a unit-gain placeholder, so its weights are proportional to photoelectrons rather than calibrated to them

This phase now asserts the opposite in three places: this docstring, the "pe" label itself, and registered_pair's "a GEDI cube SUMS to four times its stored photoelectrons". The point the comment is making — that the two sensors' weights are incommensurate and must be labelled — is exactly right and worth keeping; it just needs the same hedge the waveform notebook already uses (e.g. pe* with "proportional to photoelectrons at the store's unit-gain placeholder").

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤖 from Claude

Fixed in 7141cf4. Two of the three sites you named changed text; I checked the whole tree for a surviving contradiction and found none.

The UNITS docstring now carries the hedge explicitly:

#: What one unit of a sensor's stored WEIGHT is. ATL03 counts photons. GEDI's
#: flux is background-subtracted counts scaled by a named receiver gain, and
#: this store ships `gain: {name: "unit", value: 1.0}` -- a placeholder -- so
#: "pe" here is PROPORTIONAL to photoelectrons, not calibrated to them (the
#: published chain is ~0.4 pe/count; see `waveform_viewer.ipynb`).

registered_pair's docstring dropped the unit claim entirely — "so a GEDI cube sums to four times its stored weight", not "its stored photoelectrons".

Grepping viewers.py and both notebooks for scaled to / calibrated / photoelectrons, the only surviving assertions are the two hedged ones (viewers.py:52 above, and the waveform notebook's "this store ships a unit-gain placeholder, so its weights are proportional to photoelectrons rather than calibrated to them"), plus a generic contrast in view3d ("cells are cells, unlike photons against photoelectrons") that makes no calibration claim. The two notebooks now agree.

One residue, cosmetic and left as-is: the label itself is still bare pe rather than the pe* you suggested, so hhdc_viewer's printed lines read pe while waveform_viewer's own summary line reads pe*. No prose anywhere now claims the weights are calibrated, so this is an inconsistency in notation rather than a contradiction in fact — flagging it rather than treating it as closed.

Comment thread notebooks/hhdc_viewer.ipynb Outdated
Comment on lines +239 to +242
"\n",
"The z axis matters as much as the xy one and is easier to miss. `read_tensors`\n",
"derives its window per sensor: on this block it returns `z0 = -71.0` for ATL03\n",
"and `z0 = -59.0` for GEDI, so those two tensors are 6 bins out of register\n",

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤖 from Claude (review)

"6 bins out of register" is 12 bins for the read_tensors call this notebook actually makes.

The two z0 values are right — I get exactly -71.0 (atl03) and -59.0 (gedi) on this block. But the sentence converts the 12 m gap into bins using registered_pair's grid, not read_tensors's. view3d defaults to n_bins=256, resolution=1.0, and that is what the viewer above calls:

atl03 tensor (128, 128, 256) offset -71.0 gain 1.0
gedi  tensor  (64,  64, 256) offset -59.0 gain 1.0

So "those two tensors" are 12 bins out of register, not 6. 6 is the offset on registered_pair's own axis, which fit="degrade_resolution" pushes from the 0.5 m asked to 2.0 m (its printed line says so). registered_pair's docstring gets this right by staying in metres ("puts the two 12 m apart on this block") — the markdown could do the same, or say 12 bins at the viewer's 1 m.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤖 from Claude

Partially fixed in 7141cf4 — the number you flagged is corrected, but the same commit made the surrounding example unreproducible, so I do not think this one is closed.

The correction is right. The markdown now reads "z0 = -71.0 for ATL03 against -59.0 for GEDI, twelve bins out of register", and view3d still defaults to n_bins: int = 256, resolution: float = 1.0. I re-measured your block directly:

4331422233444 (shard 4331422233): atl03 z0=-71.0 gain=1.0 | gedi z0=-59.0 gain=1.0 | gap -12.0 m | same gain

12 m at 1 m gain is 12 bins. Your arithmetic, and the fix, both hold there.

But that block is no longer the one the notebook opens. The same commit introduced the shared shard pick (SHARD = densest_shard(...)), which now resolves to 4331422411, not 4331422233. Block 4331422233444 is in a different leaf and is never opened by a reader running the notebook top to bottom, so "on one block here" cites values that cannot be reproduced from the notebook as it stands.

On the block it actually lands on now (joint[0] = 4331422411132), and on its neighbours in the same shard:

4331422411132: atl03 z0=-201.0 gain=2.0 | gedi z0= -24.0 gain=1.0 | gap -177.0 m | GAINS DIFFER (2.0 vs 1.0)
4331422411111: atl03 z0= -34.0 gain=1.0 | gedi z0=  -9.0 gain=1.0 | gap  -25.0 m
4331422411113: atl03 z0= -23.0 gain=1.0 | gedi z0= -15.0 gain=1.0 | gap   -8.0 m
4331422411434: atl03 z0= -11.0 gain=1.0 | gedi z0=  -9.0 gain=1.0 | gap   -2.0 m

Two problems on the default block specifically: the gap is 177 m, not 12; and fit="degrade_resolution" pushes ATL03 to a 2 m gain while GEDI stays at 1 m, so "bins" is not a shared unit there at all and no single bin count is correct for both. The point the sentence is making — that read_tensors derives the window per sensor and leaves them out of register — is if anything better supported by the new block, but the concrete -71.0 / -59.0 / "twelve bins" no longer describe anything the reader will see.

Not editing this, since the fix is a judgement call about which numbers to quote: either requote against 4331422411132 (and say metres, given the gains differ), or drop the specific values and keep the qualitative point. Leaving it for the author.

@espg espg Aug 29, 2026

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤖 from Claude

Now closed in f33fecca. The partial-fix note above was right: 7141cf4c corrected the arithmetic but the same commit added densest_shard, which moved the default from shard 4331422233 to 4331422411 — so the sentence described a block the notebook no longer opens. Measured on the block it now lands on (4331422411132, at the viewer's n_bins=256, resolution=1.0):

atl03  z0 = -201.0  gain = 2
gedi   z0 =  -24.0  gain = 1

Different origin and different bin height, so the gap cannot be stated in bins at all — the earlier phrasing was unfixable as written, not just miscalculated.

Rather than substitute fresh numbers that go stale the next time the default shard moves (which is what happened here, and to the fill table on the neighbouring thread), the prose now states only the mechanism, and registered_pair prints the real windows every run:

grid  shared z = -201.0 m + bin * 4 m (DEGRADED from 0.5 m)
      — alone they would be atl03 -201.0/4 m vs gedi -24.0/0.5 m

Each sensor's solo window comes from the same chunk_z_range call applied to digests already in memory, so the contrast costs no extra I/O and cannot drift from what the code does.

One thing this surfaced that is worth recording: sharing an axis costs GEDI 8× its z resolution on this block — alone it fits the 0.5 m bins it asked for, but ATL03's deep tail drags the shared floor to -201 m and the fit policy coarsens the bin to 4 m to cover 512 m. That is inherent to demanding one axis, but it is this severe only because of an ATL03 outlier tail, so a tighter trim than the 5/95 quantiles would likely recover most of it. Left as-is and raised for review, since it changes what the export produces.

Comment thread notebooks/viewers.py Outdated
Comment on lines +210 to +213
"""Coarse sensor -> ``(fine sensor, its cells over a fine cell, its cells)``.

One GEDI o18 footprint tiles exactly 2x2 of ATL03's o19 cells, so "does
this GEDI cell see any ICESat-2 photons at all" is a block-reduce of the

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤖 from Claude (review)

The docstring asks a broader question than the code answers — and this docstring is specifically about mask semantics, so the gap reads as an error.

"does this GEDI cell see any ICESat-2 photons at all" is not what mask == 2 on 19/h_tdigest_signal reports. Per moczarr's own contract, on a signal field 1 marks the cells "whose photons were all noise" — observed ground, no signal digest. Those cells are excluded here, by design and correctly, but they were still seen by ICESat-2.

Measured on the demo block:

atl03 mask==2: 1,514 o19 cells -> 702 o18 parents; gedi over it 27 of 201 (13%)
atl03 mask>0 : 1,973 o19 cells -> 849 o18 parents; gedi over it 35 of 201 (17%)

So the printed line (and the % over atl03 pane note) says 13% where "any photons at all" would be 17%. mask == 2 is the right channel to reduce — the rest of this docstring argues that well — the question above it just needs to say signal photons / a stored digest, matching what the printed line means.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤖 from Claude

Fixed in 7141cf4. The question now names the thing mask == 2 actually reports:

One GEDI o18 footprint tiles exactly 2x2 of ATL03's o19 cells, so "does this GEDI cell sit over a stored ATL03 SIGNAL digest" is a block-reduce of the finer occupancy channel ANDed with the coarser one

"see any ICESat-2 photons at all" is gone; the rest of the docstring (the mask == 2 vs mask.astype(bool) paragraph) is unchanged and still argues the right thing.

I confirmed the distinction is still live and non-trivial on the block the notebook now defaults to (4331422411132), so the reworded question is doing real work rather than papering over a no-op:

atl03: mask==2 3,255 cells | mask>0 3,791 cells
gedi : mask==2   476 cells | mask>0   476 cells
_coverage(mask==2): {'gedi': ('atl03', 150, 476)}   -> 31.5%
_coverage(mask>0) : {'gedi': ('atl03', 162, 476)}   -> 34.0%

536 ATL03 cells here are observed-but-empty, and counting them would move the printed overlap from 150 to 162 — the same gap you measured as 27 vs 35 on the old block. The code still excludes them, and the docstring now says so.

Comment thread notebooks/hhdc_viewer.ipynb Outdated
Comment on lines +317 to +319
" print(f\"{sensor} o{order} chips \u2014 {dz:.3f} m isotropic voxels, {side}^3 = {dz * side:.0f} m cube\")\n",
" print(\n",
" f\" read {len(z):,} centroids, {wt.sum():,.0f} {UNITS[sensor]} in {read_s:.1f}s\"\n",

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤖 from Claude (review)

Minor, and it is the one place this phase's own honest-reporting rule slips: the read line prints wt.sum() — the total read, not the total written. At the default o22 they are the same, but at the o24 the cell right below invites:

atl03 o24 chips — 0.389 m isotropic voxels, 128^3 = 50 m cube
  read   22,252 centroids, 22,252 ph in 14.4s  (2,565 above the z window)
  wrote  298 of 1024 chips (726 skipped empty), 16,717 voxels filled, in 3.5s

11.5% of the photons are not in the file, and the only number that says so is a centroid count in parentheses on the read line. registered_pair sets the better precedent — its totals come off cube.sum(), so they are what is on disk by construction. Adding the binned weight to the wrote line (e.g. 19,687 of 22,252 ph binned) would close it; the drop accounting itself is correct — I confirmed z0 = floor(chip min) makes negative iz impossible, that a written chip always holds at least one voxel, and that nothing else is dropped silently.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤖 from Claude

Fixed in 7141cf4, and I checked the arithmetic rather than the wording. The wrote line now carries the binned weight, with the drop stated as a difference:

print(f"  wrote  {len(meta)} of {n} chips, {kept:,} of {wt.sum():,.0f} {UNITS[sensor]} "
      f"({wt.sum() - kept:,.0f} outside a chip's z window), ...")

where kept accumulates chip.sum() — post-rint, post-clip, so it is what is on disk by construction, the registered_pair precedent you pointed at.

Live output on the notebook's current default block (4331422411132):

o22:  read 89,205 centroids, 89,205 ph | wrote 46 of 64 chips, 87,218 of 89,205 ph (1,987 outside a chip's z window)
o24:  read 89,205 centroids, 89,205 ph | wrote 480 of 1024 chips, 87,833 of 89,205 ph (1,372 outside a chip's z window)

I then reimplemented the binning independently and counted the drops by cause, to check the parenthetical is neither double-counting nor hiding a second drop path:

all weights exactly 1.0? True  (min=1.0 max=1.0)  -> np.rint loses nothing
o22: photons assigned to some tile = 89,205 (== read: True)
     dropped by iz>=n_bins = 1,987   negative iz = 0
     read - kept = 1,987   MATCH: True
o24: photons assigned to some tile = 89,205 (== read: True)
     dropped by iz>=n_bins = 1,372   negative iz = 0
     read - kept = 1,372   MATCH: True

So read - written equals the iz >= n_bins weight exactly at both orders: every photon lands in some tile (skipped-empty chips drop nothing, since a skipped tile had no photons), z0 = floor(min) keeps iz non-negative as you'd confirmed, and rounding contributes zero because every ATL03 centroid weighs 1.0. The label "outside a chip's z window" is the whole of the difference, not part of it.

One scope note: the rint-is-free half of that holds because this store's ATL03 weights are unmerged. voxel_chips(sensor="gedi") would have fractional weights, where read - written would fold rounding in with the z clip and the parenthetical would be slightly over-attributed. The notebook only ever calls it for ATL03, so nothing printed is wrong today.

@espg espg left a comment

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤖 from Claude (review)

Adversarial review of f33fecca..e59ce041 (export.py + fit_window, viewer_stats, paired_blocks, the stats= signature change, and the export prose). Six findings, all reproduced — two on the public store's own block 4331422411132. viewer_stats and paired_blocks are behaviourally identical to the notebook code they replaced; nothing to raise there.

Comment thread notebooks/export.py Outdated
zs, ws = np.asarray(z)[order], np.asarray(wt)[order]
span = n_bins * dz
if zs[-1] - zs[0] <= span: # fits as it stands
return float(np.floor(zs[0])), 0.0, 0.0, 1.0

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤖 from Claude (review)

The fast path drops the top of the data while reporting kept=1.0, trim_low/high=0.0.

z0 = float(np.floor(zs[0])) can sit up to 1.0 m below zs[0], but the window voxel_chips then bins into is [z0, z0 + span). So the guard zs[-1] - zs[0] <= span is not the condition that matters — zs[-1] < z0 + span is. Everything between z0 + span and zs[-1] fails the iz < n_bins mask at export.py:118 and is dropped, while meta.json records kept: 1.0, trim_low: 0.0, trim_high: 0.0 and the worst kept line (which filters mm["kept"] < 1.0) stays silent about it.

The loss is 1 m / dz bins, so it scales with order. Synthetic, 100k uniform points spanning 99.9% of the window starting at z=100.9:

order   dz (m)   span (m)   reported kept   ACTUAL kept   bins lost
   22  1.55440    198.963           1.000        0.9965           1
   24  0.38860     49.741           1.000        0.9829           3
   26  0.09715     12.435           1.000        0.9285          10
   29  0.01214      1.554           1.000        0.4214          74

Confirmed on the store, block 4331422411132, 19/h_tdigest_signal (89,205 photons), replaying the exact per-chip loop:

o22:    46 chips |    0 chips report kept=1.0 while dropping     0 ph
o24:   480 chips |    2 chips report kept=1.0 while dropping     2 ph
o26:  3255 chips |   33 chips report kept=1.0 while dropping    38 ph
o29: 25584 chips | 1767 chips report kept=1.0 while dropping 3,172 ph

The default o22 run is clean, but the notebook's own commented-out order=24 line already trips it, and the markdown two cells up advertises "down to 1.21 cm at order 29" — where reported kept is 0.727 against an actual 0.692. This is the exact thing the markdown promises does not happen ("nothing is clipped silently") and the docstring calls "worse but never silent".

Cheapest fix is to make the fast path honour the window it actually hands back:

z0 = float(np.floor(zs[0]))
if zs[-1] < z0 + span:  # fits as it stands
    return z0, 0.0, 0.0, 1.0

leaving the trim path (whose kept is computed from the data and is honest at the top) to handle the residual case.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤖 from Claude

Fixed in a076d578 — verified independently, including the o29 case.

The guard now computes the floor first and tests the window that is actually handed back:

total = float(ws.sum())
z0 = float(np.floor(zs[0]))
if total <= 0 or zs[-1] < z0 + span:
    return z0, 0.0, 0.0, 1.0

Your synthetic, replayed against HEAD (100k uniform points spanning 99.9% of the window, base z=100.9). Reported kept now equals what voxel_chips actually bins at every order — the fast path no longer fires, and the trim path's kept is honest:

order   dz (m)   span (m)   reported   ACTUAL   was reported
   22  1.55440    198.963     0.9964   0.9964         1.000
   24  0.38860     49.741     0.9829   0.9829         1.000
   26  0.09715     12.435     0.9287   0.9287         1.000
   29  0.01214      1.554     0.4201   0.4201         1.000

On the store, block 4331422411132, 19/h_tdigest_signal (89,205 photons — your number), replaying the exact per-chip loop:

order  chips  divergent chips  ph lost silently   reported kept   actual kept
   22     46                0                 0           0.999         0.999
   24    480                0                 0           0.994         0.994
   26   3255                0                 0           0.893         0.893
   29  25584                0                 0           0.690         0.690

The o29 case you measured — 1,767 chips reporting kept=1.0 while dropping 3,172 ph, aggregate 0.727 reported against 0.692 actual — is gone: zero divergent chips, zero photons lost silently, and reported now tracks actual at 0.690. (0.690 rather than your 0.692 because the _wq midpoint fix in the same commit shifts the trim quantiles by half a centroid, so the fitted windows differ slightly.)

The guard fix is load-bearing on its own. Fuzzing 12,000 synthetic chips across orders 22/24/26/29 and five shapes (uniform overflow, ground cluster + sparse canopy, bimodal give-up, all-equal z, heavy low tail):

HEAD (new guard + np.floor)      :    0 chips where reported kept != actual
revert the guard only            :   95 chips diverge
revert the np.floor only         : 1505 chips diverge

So neither fix is redundant — each catches cases the other does not.

Edge cases also check out (single centroid, all-equal z, repeated ties with span overflow, a chip whose low tail sits entirely under the floor): reported equals actual in all of them, and the worst kept diagnostic now fires because kept < 1.0 is reached honestly.

Comment thread notebooks/export.py Outdated
if not len(m):
continue
z0, lo_d, hi_d, frac = fit_window(z[m], wt[m], n_bins, dz)
iz = ((z[m] - z0) / dz).astype(np.int64)

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤖 from Claude (review)

.astype(np.int64) truncates toward zero, so data BELOW the fitted floor is kept and dumped into bin 0.

For z0 - dz < z < z0 the quotient is in (-1, 0) and truncation gives iz == 0, which passes iz >= 0 on the next line. Those centroids are then accumulated at acc[..., 0] — kept, at the wrong height.

The pre-refactor version was safe by construction: z0 = np.floor(z[m].min()) guaranteed iz >= 0, and the mask was only iz < n_bins. fit_window removed that guarantee — when it trims a low tail, z0 sits above the data minimum on purpose — but the mask picked up a >= 0 test rather than a floor, so the sub-floor slab is silently folded in instead of being dropped.

Synthetic (100k centroids, a tight low cluster below the 2% cut plus canopy above):

z0 = 408.000  lo_drop=0.07  hi_drop=0.05  reported kept = 0.93074
actual kept by voxel_chips        = 0.93494
below-floor points kept           = 420 (0.42% of the chip), all forced into bin 0
  their true z spans [406.454, 407.998] -- up to 1.546 m below the floor
np.floor would have given iz = -1 for them

Confirmed on the store, block 4331422411132: 9 ph at o22, 29 ph at o24, 90 ph at o26 land in bin 0 from below the floor, and reported kept and actual kept diverge in the same direction at every order.

The bound on the contamination is "everything in [z0 - dz, z0)" — one full voxel of relief, wherever the low tail happens to be dense right under the cut.

iz = np.floor((z[m] - z0) / dz).astype(np.int64)

Secondary, same lines: trimmed += int(wt[m][~keep].sum()) truncates a float sum per chip, and np.rint at line 123 rounds per occupied voxel, so written + trimmed does not close against read even once the mask is fixed. The print at line 145 invites exactly that subtraction.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤖 from Claude

Primary point fixed in a076d578; the secondary accounting point is NOT fixed and still stands.

The cast now floors first, with the reasoning kept in a comment at the site:

iz = np.floor((z[m] - z0) / dz).astype(np.int64)

Sub-floor centroids are now excluded, not mis-binned. Over 12,000 fuzzed chips (orders 22/24/26/29, five shapes) I found 6,248 chips holding 389,703 centroids below z0. At HEAD I asserted per chip that none of them survive the iz >= 0 mask — none did. The old truncating cast would have folded 3,387 of those into bin 0 at the wrong height.

On the store, block 4331422411132, the photons the old cast would have dumped into bin 0 from below the floor: 9 at o22, 42 at o24, 108 at o26, 56 at o29 — your 9/29/90 shifted slightly because the _wq fix in the same commit moves the fitted windows. At HEAD, reported kept and actual kept agree to 1e-9 on every one of the 46/480/3,255/25,584 chips, which is the same statement: fit_window's kept counts (zs >= z0) & (zs < z0 + span) and the floored mask keeps exactly that set, so the two can only agree if nothing sub-floor is kept.

The floor fix is load-bearing on its own — reverting it while keeping the new fast-path guard reintroduces divergence on 1,505 of the 12,000 fuzzed chips (reverting the guard alone: 95). Neither is redundant.

Still open — the closure nit in your last paragraph. a076d578 did not touch trimmed += int(wt[m][~keep].sum()) or the np.rint at what is now line 140, so written + trimmed still need not close against read. It happens not to bite on the demo block, because ATL03 signal-photon weights there are integral — it closes exactly at all four orders:

read = 89,205 photons
order    written   trimmed        w+t   read-(w+t)
   22     89,154        51     89,205            0
   24     88,702       503     89,205            0
   26     79,675     9,530     89,205            0
   29     61,521    27,684     89,205            0

But it is latent for fractional weights, which is what the GEDI rx_flux path carries. One synthetic chip, 20k centroids with weights drawn from U(0.01, 3.0):

read = 30,269.80   written = 19,991   trimmed = 10,255   w+t = 30,246   gap = 23.80  (0.08%)

Leaving that standing rather than folding it: it is a reporting-precision issue in a print, not data loss, and the fix touches how trimmed and written are accumulated — a judgement call about whether the printed line should carry fractional weight. Flagging it for the author rather than deciding it here.

Comment thread notebooks/export.py Outdated

z0 = float(np.floor(med - span / 2)) # give up trimming; centre on the mass
inside = float(ws[(zs >= z0) & (zs < z0 + span)].sum() / cum[-1])
return z0, np.nan, np.nan, inside

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤖 from Claude (review)

The give-up branch returns np.nan, which json.dumps writes into meta.json as a bare NaN token — not valid JSON.

lo_d/hi_d go straight into meta[name]["trim_low"]/["trim_high"] at export.py:132-133, and meta.json is the manifest shipped inside the exported .npz:

>>> json.dumps({"chip": {"z0": -100.0, "dz": 1.5543972149675078,
...                      "trim_low": nan, "trim_high": nan, "kept": 0.5}})
'{"chip": {"z0": -100.0, "dz": 1.5543972149675078, "trim_low": NaN, "trim_high": NaN, "kept": 0.5}}'

Python's json.loads accepts it by default, but JSON.parse, jq, serde_json and json.loads(..., parse_constant=...) all reject a bare NaN. Reproduced with a bimodal chip (two clusters 5,000 m apart, neither trimmable inside max_drop), which is exactly the kind of chip this branch exists for.

None (→ null) reads the same to a consumer and stays parseable; the realized drop fractions would be better still, since the docstring above promises "the weight fraction trimmed off each end" and NaN is not one.

Related, same expression: a chip whose weights are all zero makes cum[-1] == 0, and line 78 raises RuntimeWarning: invalid value encountered in scalar divide and returns kept = nan as well.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤖 from Claude

Fixed in a076d578 — both halves, verified.

The give-up branch returns None, and the zero-total case is guarded before any division:

total = float(ws.sum())
...
if total <= 0 or zs[-1] < z0 + span:
    return z0, 0.0, 0.0, 1.0
...
return z0, None, None, inside  # None, not NaN: `meta.json` has to parse

_wq grew the matching guard (if cum[-1] <= 0: return float(zs[0])).

Give-up branch round-trip. Your bimodal chip (two clusters 5,000 m apart, neither trimmable inside max_drop) reaches the branch and produces:

{"chip": {"z0": -97.0, "dz": 1.5543972149675078, "written": 123,
          "trim_low": null, "trim_high": null, "kept": 0.5018663236870695}}
  • json.loads(..., parse_constant=<raises>) — passes (it would raise on a bare NaN).
  • JSON.parse (node) — parses. Confirmed the same call rejects {"a":NaN} with Unexpected token 'N'.
  • jq — parses. One correction to the finding: jq 1.7.1 is not strict here — it accepts a bare NaN and silently coerces it to null, so jq would not have caught the old output. JSON.parse and a parse_constant hook do; I did not test serde_json.

Zero-weight chip. Ran fit_window on 50 centroids with all-zero weights under warnings.simplefilter("error"), so a RuntimeWarning would have been raised as an exception. None was. It returns z0=0.0, trim=(0.0, 0.0), kept=1.0 — a finite kept, and json.dumps of that meta contains no NaN. _wq on zero weights returns zs[0] rather than dividing.

Sanity-checked the same across the fuzz set and the edge cases (single centroid, all-equal z, ties, entire low tail under the floor) — every meta.json I generated serialised without a NaN or Infinity token.

Comment thread notebooks/export.py Outdated

def _wq(zs, cum, q):
"""Weighted quantile of pre-sorted `zs` given the cumulative weight `cum`."""
return float(np.interp(q * cum[-1], cum, zs))

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤖 from Claude (review)

_wq interpolates on the raw np.cumsum, so every quantile is biased low by half a centroid's weight.

The standard weighted-quantile plotting position is (cum - ws/2) / cum[-1]; using cum directly puts the whole of zs[i]'s weight below zs[i]. On n equally-weighted centroids uniformly spanning 0–100 m:

n=  2  _wq median =  0.000   true  50.000   bias -50.000 m
n=  3  _wq median = 25.000   true  50.000   bias -25.000 m
n=  5  _wq median = 37.500   true  50.000   bias -12.500 m
n= 11  _wq median = 45.000   true  50.000   bias  -5.000 m
n=101  _wq median = 49.500   true  50.000   bias  -0.500 m

Two places it bites, both in this file:

  1. The give-up branch centres the fixed window on med (line 77). A sparse chip — and at o24/o26 most chips are sparse — gets the window mis-centred by a large fraction of its own relief, which is the one branch that has no other defence.
  2. np.interp left-clamps below cum[0], so when the bottom centroid holds more weight than lo_d, the low cut moves lo not at all while meta.json still records trim_low: 0.02. Reproduced with a bottom centroid holding 4% of the weight at z=50 against a mass at 400–900 m: _wq(q) returns 50.0 for every q <= 0.04.

np.interp(q * cum[-1], cum - ws / 2, zs) fixes both.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤖 from Claude

Fixed in a076d578 — verified.

_wq now takes the weights rather than a precomputed cumsum, and interpolates on midpoint positions:

cum = np.cumsum(ws)
if cum[-1] <= 0:
    return float(zs[0])
return float(np.interp(q * cum[-1], cum - 0.5 * ws, zs))

Your bias table, rerun against HEAD — n equally-weighted centroids uniformly spanning 0–100 m:

n=  2  _wq median =  50.000   true 50.000   bias 0.000 m   (was -50.000)
n=  3  _wq median =  50.000   true 50.000   bias 0.000 m   (was -25.000)
n=  5  _wq median =  50.000   true 50.000   bias 0.000 m   (was -12.500)
n= 11  _wq median =  50.000   true 50.000   bias 0.000 m   (was  -5.000)
n=101  _wq median =  50.000   true 50.000   bias 0.000 m   (was  -0.500)

Exact at every n, so point (1) — the give-up branch mis-centring a sparse chip on med — is resolved at the source.

Point (2), the left clamp. Your case, a bottom centroid holding 4% of the weight at z=50 against a mass at 400–900 m:

_wq(q=0.00) =  50.000
_wq(q=0.01) =  50.000
_wq(q=0.02) =  50.000
_wq(q=0.04) = 330.000   (was 50.000 -- the cut now moves)
_wq(q=0.05) = 402.632
_wq(q=0.10) = 428.947

The 2% first cut is the residual clamp, and I read that as correct rather than a leftover bug: with midpoint positions the bottom centroid's own position sits at w[0]/2, i.e. q=0.02, so any q at or below that is genuinely below the first plotting position and clamping to zs[0] is what the convention prescribes. The clamped interval is now exactly half a centroid wide instead of a whole one, which is the improvement the finding asked for.

Downstream effect worth noting: because the trim quantiles moved, the fitted windows on the store shifted slightly — the aggregate o29 kept on block 4331422411132 is 0.690 at HEAD where the finding measured 0.692 actual against the old _wq.

Comment thread notebooks/hhdc_viewer.ipynb Outdated
"Both are the user's call, i.e. client-side functions. Because the original locations\n",
"are stored, a tensor can be built on a 12 m, 1.5 m or 0.38 m spatial grid — down to\n",
"1.21 cm at order 29. In the examples below, range bins are isotropic with the\n",
"spatial grid by default, or any resolution you ask for when exporting jointly.\n",

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤖 from Claude (review)

"or any resolution you ask for when exporting jointly" is falsified by the very cell this paragraph introduces.

registered_pair hardcodes fit="degrade_resolution" (export.py:181), which is the moczarr policy that doubles resolution until the trimmed range fits n_bins — the requested resolution is a floor, not a promise. I ran the notebook's own default call against the public stores, block 4331422411132:

registered o19 — (128, 128, 128) float32 each, 54.1s read
  grid   shared z = -201.0 m + bin * 4 m (DEGRADED from 0.5 m) — alone they would be
         atl03 -201.0/4 m vs gedi -24.0/0.5 m

Asked for 0.5 m, got 4 m — 8x coarser — on the demo's default block. The code is fine and even prints DEGRADED from 0.5 m; it is the prose that overclaims, and a reader who takes the sentence at face value will be surprised by the cell's own output three lines later.

Something like "…or a resolution you name when exporting jointly, coarsened by powers of two if the pair's relief will not fit" would match what the cell prints.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤖 from Claude

Fixed in a076d578 — the prose no longer promises the resolution.

Cell 8 of hhdc_viewer.ipynb at HEAD now reads:

In the examples below, range bins are isotropic with the spatial grid by default; the joint export takes whatever z resolution you ask for, but coarsens it when the two sensors' shared relief will not fit — and prints that it did.

replacing "…by default, or any resolution you ask for when exporting jointly."

Verified the second clause is true of the code rather than a second unbacked claim — registered_pair does print it, at export.py:229:

fit = "as asked" if abs(dz - resolution) < 1e-9 else f"DEGRADED from {resolution:g} m"

so the demo's default call, which the finding measured coming back at 4 m against a requested 0.5 m, prints DEGRADED from 0.5 m and the paragraph now sets that expectation instead of contradicting it. I confirmed the text and the print statement at HEAD; I did not re-run the paired export against the stores, so the 8x figure itself is still your measurement rather than mine.

fit="degrade_resolution" at export.py:198 is unchanged, which is right — the finding was scoped to the prose, not the policy.

One incidental note on the same commit: it rewrote every em-dash in this notebook from a literal to a escape. Content-identical, but it makes the diff on this file much noisier than the one-sentence change it carries.

Comment thread notebooks/viewers.py Outdated
found = {w for per in tally.values() for w in per}
if not found & within:
raise ValueError(
f"the tally holds no block beneath shard {shard} — `handles`/`tally` and "

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤖 from Claude (review)

Stale after 7ebed910: the message names tally, but view3d's third argument is now stats. Someone who hits this goes looking for a tally= keyword that no longer exists — `handles`/`stats` would point them at the argument they actually passed.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤖 from Claude

Fixed in a076d578 — verified.

viewers.py:509 at HEAD:

raise ValueError(
    f"the stats hold no block beneath shard {shard} — `handles`/`stats` and "
    f"`shard` disagree. Open and view the SAME shard (one name, used twice)."
)

Checked it against the signature it is describing: view3d takes stats=None (viewers.py:392) and unpacks it at :429–430 (tally = (stats or {}).get("tally")), so stats is the argument a caller actually passes and the message now names it. tally survives only as a local for the inner dict — which is also the key viewer_stats returns at :172 ({"tally": ..., "joint": ...}), so the name is not stale in those places, just no longer part of the caller-facing surface. Grepped the rest of the module for other user-visible tally references in messages or docstrings: the remaining hits are all locals or that dict key, none in an error string.

@espg
espg marked this pull request as ready for review August 29, 2026 21:39
@espg
espg merged commit 5800986 into main Aug 29, 2026
6 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant