Repository navigation
Conversation
subset.bounding_box asked for Grid.bounds, which builds a per-face box for every face in the mesh -- materializing the node coordinate and connectivity arrays to do it -- and only then tested which faces fall inside the region. The cost was proportional to the mesh rather than the crop, so cropping a few hundred faces out of a large grid needed memory proportional to the whole grid. On a 300M-face SCRIP mesh that exceeds the machine, and the subset is not slow but impossible (#1778). Face latitudes are already resident, so they now discard faces that cannot be in the result before any bounds are computed; exact bounds are built for the survivors only. Measured on a 12-corner SCRIP mesh, crop ~0.007% of the faces: 1,600,000 faces 0.84 GiB, 2.6 s -> 0.00 GiB, 0.8 s 3,200,000 faces 1.48 GiB, 5.1 s -> 0.00 GiB, 1.6 s Two things measurement settled that reasoning got wrong. The first attempt filtered longitude on face centres too, with a margin sized from the mean face area; that is unsafe near the poles, where a face is recorded as spanning every longitude and its bounds sit up to 354 degrees from its own centre on a HEALPix z3 grid. No finite margin covers that, so longitude is left to the exact stage. The second attempt kept a latitude margin for the same reason -- and it is unnecessary, because faces_within_lat_bounds keeps a face only when its bounds are fully contained in the interval, despite the docstring saying "overlap". A contained face has a contained centre, so the filter is conservative by equality rather than by a fudge factor. The exact path is unchanged below a million faces and whenever bounds are already computed, since the filter cannot repay itself there. Tests compare the filtered and exact paths across three resolutions and eight boxes including antimeridian-spanning and polar cases, and pin the two measurements the design rests on: that longitude bounds reach far beyond a face centre, and that the latitude test is containment rather than overlap. If either stops being true the pre-filter is wrong and these fail.
The latitude prefilter alone leaves a global belt: every face at the crop's latitude, at every longitude. On a HEALPix mesh a Texas-sized box still kept 8.3% of the faces. Longitude needs a margin where latitude did not, and the margin has to widen towards the poles because meridians converge -- a fixed-area face spans more longitude the higher it sits. Normalising the measured reach by 1/cos(lat) collapses it to a constant 4.39 face widths across zooms 4, 5 and 6, which is what makes the bound derivable rather than tuned. Eight face widths, latitude-corrected, had zero violations from zoom 3 through 7. Two kinds of face cannot be filtered on a centre at all and are always kept: one whose bounds are stored min > max crosses the antimeridian, so its centre is not between them, and one containing a pole is recorded as spanning every longitude. Both are under 1% of faces by zoom 6. Texas candidates fall from 8.3% to 0.5% of the mesh, improving with resolution. 44 of 44 boxes return exactly what the unfiltered path returns, including antimeridian-spanning and polar queries.
Sevans711
added
scalability
ASV BenchmarkingBenchmark Comparison ResultsBenchmarks that have improved:
Benchmarks that have stayed the same:
Benchmarks that have got worse:
|
Sevans711
requested changes
Sep 23, 2026
Collaborator
There was a problem hiding this comment.
This is a great idea. It definitely still needs some work before it is ready to merge; see inline comments. Two more notes here:
(1) This doesn't fully close the original issue; #1778 is scoped to include the face_node_connectivity issue too. The original message in this PR should be updated accordingly, to avoid the "Closes #XXX" pattern. (EDIT: actually, upon a more careful read through of 1778, it's a bit unclear to me what the intent was. The title seems to indicate it is scoped just to the bounding-box subset issue, but the text seems to include the connectivity. If you meant to scope it only to the former, please clarify, and then feel free to keep this PR marked as "would close 1778"!)
(2) The ASV benchmarks suite does not show any relevant improvements. Can benchmarks be added which show improvement from this PR? Ideally, at least one showing memory improvement, and at least one showing speed improvement. That would help prevent performance degradation in the future.
Each edge contributed either its interior latitude extreme or its first node, never both, so a face whose edges all bulge poleward lost its lowest node from its own bounds. geoflow-small has 70 such faces.
A face's bounds contain each of its nodes, so a face can only lie inside the box if all its nodes do. Bounds are then computed for those candidates alone, on the nodes they use, so the cost follows the region rather than the mesh. Replaces the centre-based filter and its tuned margins.
rajeeja
changed the title
Contributor
Author
|
@Sevans711 Reworked in 7e577f7..dd351b0: the centre filter is replaced by a node-coordinate screen with no margins, the bounds kernel now includes every node (it was dropping some), tests rewritten to check what they say, and |
Sevans711
requested changes
Sep 28, 2026
Sevans711
left a comment
Sevans711
left a comment
Collaborator
There was a problem hiding this comment.
Getting much closer! Some small inline comments.
The ASV benchmarks look good and show the relevant changes, which addresses my previous comment (2).
For my previous comment (1), your updated description is actually even more confusing to me. Do you think the connectivity point is in scope or out of scope for the original issue? (in scope --> this PR shouldn't close the issue. out of scope --> this PR should close the issue.)
| def test_face_bounds_contain_every_node(gridpath): | ||
| """A face's latitude bounds must not be narrower than its own nodes. | ||
|
|
||
| Each edge used to contribute either its interior latitude extreme or its |
Collaborator
There was a problem hiding this comment.
What do you mean by "used to"? Is this relative to main, or relative to an earlier commit in this PR?
Contributor
Author
There was a problem hiding this comment.
Relative to main — the bounds kernel on main has this bug today. Each edge contributed either its interior latitude extreme or its first node, never both, so a face whose edges all bulge poleward lost its lowest node (70 faces in geoflow-small). Fixed in this PR's first commit, 7e577f7.
It is in scope here because the screen relies on a face's bounds containing its nodes; without the fix the screen could drop a face the exact test keeps. Docstring now says all of that explicitly.
Collaborator
There was a problem hiding this comment.
Ah, I see. My concern here was mostly just that after this PR merges, language like "on main" and "used to" (referring to this PR versus main) and "this PR" won't really make sense anymore, because the test will become part of main. Could you update the language in the docstring so that it will still make sense after merging? This could be as simple as swapping those phrases to something like "before PR 1779" and "after PR 1779" or "fixed in PR 1779" instead.
…ian logic - Move bounds/intersections/coordinates imports to module scope; no circular import. Monkeypatch targets in the tests move with them. - Add a docstring to each test in TestBoundingBoxSubset saying what it pins and why. - Spell out the lon_min > lon_max case in _faces_with_nodes_within_box: the box is the union of two intervals, so a node is outside only in the gap between them. - test_face_bounds_contain_every_node: say the behaviour it pins is relative to main, fixed in this PR's first commit.
Contributor
Author
|
Updated with latest Locally: @Sevans711 ready for another look when you have time. |
Sevans711
reviewed
Oct 20, 2026
Sevans711
left a comment
Sevans711
left a comment
Collaborator
There was a problem hiding this comment.
Looking great overall. Still one minor nit (see the last remaining unresolved inline comment from prior review).
My only other question now is whether the increase to the nogil_scaling.GILScaling.track_gil_scaling benchmark is something we should be concerned about. @cmdupuis3 do you have any insights about what to expect from this benchmark in this PR? Is it an issue with the PR itself, or more likely to be noise?
Aside from those points, I think I would be ready to approve!
cmdupuis3
requested changes
Oct 28, 2026
cmdupuis3
left a comment
cmdupuis3
left a comment
Collaborator
There was a problem hiding this comment.
Overall it seems good, the improvements are verified. I think there is more optimization work that can be done, but I think there are a few tradeoffs of peakmem vs speed that I didn't want to impose on this PR. Basically, np.intersect1d dominates in _faces_in_bounding_box, and there's a related np.unique that implements a sort we don't need. I'm also curious if some of the caching on main can be salvaged if the bounding box is relatively large compared to the grid.
| return bounds_lon | ||
|
|
||
|
|
||
| @njit(cache=True, parallel=True, nogil=True) |
Collaborator
There was a problem hiding this comment.
Suggested change
| @njit(cache=True, parallel=True, nogil=True) | |
| @parallel_njit |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Closes #1778 (its first point; the
face_node_connectivitycost is a reader problem and stays open in #1778's discussion).Overview
subset.bounding_boxcomputedGrid.boundsfor every face before testing which fall inside the box, so a small crop of a large mesh cost memory and time proportional to the mesh.It now screens faces on their node coordinates first: a face's bounds contain each of its nodes, so a face can only lie inside the box if all its nodes do. Bounds are computed for those candidates only, on the nodes they use. The screen needs no margin and makes no assumption about face size or shape. The result is identical to the whole-mesh path, and cached bounds are used when present.
Along the way, the bounds kernel dropped a node from its own face's bounds when both adjacent edges bulge poleward (70 faces in
geoflow-small); fixed in the first commit, since the screen relies on bounds containing the nodes.faces_within_lon_bounds/faces_within_lat_boundsdocstrings now say what they do (containment, not overlap).Measured
12-corner SCRIP mesh from the #1778 MCVE, crop of ~0.007% of faces:
HEALPix zoom 8 (786k faces), Texas box: 88 MB / 973 ms → 13 MB / 11 ms.
Testing
TestBoundingBoxSubsetcompares the screened path against the whole-mesh path on 7 grids × 12 boxes (polar, antimeridian-crossing, whole globe, empty), assertsboundsis never populated, counts the faces the bounds kernel sees, checks cached bounds are used, that faces straddling a box edge are excluded, nodes on the box edge, antimeridian and pole faces explicitly, and an empty region on a partial grid.test_face_bounds_contain_every_nodepins the kernel fix.New
benchmarks/subset.py:time_andtrack_peakmem_forbounding_boxon fresh oQU and HEALPix grids.PR Checklist
General
Testing & Benchmarking
benchmarks/subset.py)Documentation and Examples
docs/api.rstAI Disclosure
AI Usage: Claude