Better chunking strategies for constlat intersections and zonal routines. - #1624
Open
cmdupuis3 wants to merge 76 commits into
Open
Better chunking strategies for constlat intersections and zonal routines.#1624cmdupuis3 wants to merge 76 commits into
cmdupuis3 wants to merge 76 commits into
Conversation
…rite intersections, add 241 baseline testsgit status! - most came from accusphere
…us/mask, dispatcher
- benchmarks/geometry_kernels.py: ASV micro-benchmarks for all three layers of the EFT intersection stack (_accux_gca, _try_gca_gca_intersection, gca_gca_intersection, _accux_constlat, _try_gca_const_lat_intersection, gca_const_lat_intersection) plus EFT primitives and point-in-polygon; all functions warmed before timing so results reflect steady-state cost - test/test_plot.py: add test_to_raster_auto_extent verifying that the axis limits change and the raster contains finite data
…rectness fixes Review comments addressed: - Remove "near-double precision" / "sufficient" overclaims; say "roughly twice as accurate" and note the robustness tier boundary clearly - Explain _lon_bounds_from_vertices is required for UXarray antimeridian encoding and cannot be removed - Add block comment before _no_extreme functions clarifying they are pre-existing edge screeners unrelated to the EFT stack - Document SoS as explicit future work in _point_in_polygon_sphere docstring - L2 pos_fin/neg_fin: replace ternary with int(); exploit neg=-pos symmetry - Label computation: drop dead local*0 term, use integer mask arithmetic - Remove vertex-lat snap from bounds: _face_location_info already captures interior arc extrema accurately via the compensated kernel - _ON_MINOR_ARC_TOL: document intentional 1e-10 vs C++ 1e-8 divergence Bug fixes: - on_minor_arc: add antipodal-endpoint guard; a x b = 0 for antipodal inputs so every point on the great circle passes the collinearity test (false pos) - bounds.py: replace mask arithmetic use_ext*z_ext + (1-use_ext)*z_edge with plain if/else; 0*NaN = NaN propagates when norm=0, if/else does not - _point_in_polygon_sphere: ray-nudge now restarts the loop from i=0 so all edges are counted with the same ray (mid-loop nudge corrupted crossing parity) Cleanup: - Remove _flip_sign, _SIGN_NEG, _SIGN_POS, _SIGN_ZERO dead code from point_in_face.py; inline literals in _counts_as_crossing - Remove _SNAP_TOL_DEG constant and snap_tol_deg parameter throughout bounds.py - Notebook: fix Grid.get_point_on_face -> get_faces_containing_point; remove incorrect geometry.py row from Section 4 table; add accucross_pair and acc_sqrt_re to Section 2 building-blocks table
…PI name - ci/environment.yml: pin tornado<6.5.7 to avoid ssl.SSLError in panel 1.9.3 on Python 3.11 Windows (conda-forge regression, 2026-06-10) - intersections.py: remove _gca_gca_intersection_cartesian shim (dead code); add comment explaining _snap_const_lat_endpoint snap_sq constant - test_intersections.py: update 4 call sites to use gca_gca_intersection directly - spherical-geometry-accuracy.ipynb: fix stale Grid.get_point_on_face -> Grid.get_faces_containing_point (2 occurrences)
…rite intersections, add 241 baseline testsgit status! - most came from accusphere
…us/mask, dispatcher
Reconcile diverged accusphere branch. Resolutions: - intersections.py: restore inline=always on L1 kernels (_accux_constlat, _accux_gca) for allocation scalar-replacement - point_in_face.py: keep restart-loop ray casting (consistent parity), adopt named sign constants, drop unused _flip_sign - arcs.py: keep antipodal-endpoint guard in on_minor_arc - bounds.py: keep vertex-latitude snapping (snap_tol_deg) path - computing.py: keep detailed docstring with SIAM/EGUsphere references
Add an LLVM fma intrinsic and route two_prod through a single fused multiply-add for its error term on hardware that supports it, selected at import time and validated to be bit-exact against the Veltkamp split. Falls back to the portable Veltkamp form otherwise, so there is no hard FMA dependency. The FMA path is ~2x faster in the compensated geometry kernels (each two_prod drops from ~17 flops to one FMADD) and is numerically identical: all 241 AccuSphGeom baseline cases pass unchanged.
Add _accux_constlat_scalar, which takes the arc endpoints as six scalars and returns the candidate coordinates as scalars instead of two np.empty(3) arrays. _accux_constlat now wraps it so the array API is unchanged. Returning scalars lets Numba keep the candidates in registers, so a batch loop over many edges does no per-point heap allocation. On a 16M-point const-lat sweep this is ~2.7x faster than the array-returning path and drops the AccuX/FP64 cost ratio from ~19x to ~7x. Bit-identical results; all 241 AccuSphGeom baseline cases pass.
…_lat_intersection
Collaborator
Author
|
pre-commit.ci autofix |
for more information, see https://pre-commit.ci
cmdupuis3
added
the
benchmarking
ASV BenchmarkingBenchmark Comparison ResultsBenchmarks that have improved:
Benchmarks that have stayed the same:
Benchmarks that have got worse:
|
cmdupuis3
marked this pull request as ready for review
July 29, 2026 22:38
Contributor
|
Thanks for putting this together. I'll take a look at this and post my comments soon |
rajeeja
requested changes
Jul 30, 2026
rajeeja
left a comment
rajeeja
left a comment
Contributor
There was a problem hiding this comment.
Verified independently: subset builder is bit-identical to whole-grid builder (tested HEALPix + mixed-polygon MPAS), vectorized screeners match brute-force reference, both dask/AttributeError bugs reproduce on main and are fixed here. Two small comments below, non-blocking on correctness but worth addressing before merge.
| nedge = n_nodes_per_face[f] | ||
| weights[pos] = _compute_band_overlap_area( | ||
| faces_edge_nodes_xyz[f, :nedge], zmin, zmax | ||
| if partial.size: |
Contributor
There was a problem hiding this comment.
New branch (zero partial faces). Add a test for a band with only fully-contained faces.
Collaborator
Author
There was a problem hiding this comment.
Not sure what you mean here... The weights are already computed for all faces overlapping the band, this branch only reweights partial faces, so a strict fully-contained set of faces is already accounted for by the time this branch is reached.
Contributor
There was a problem hiding this comment.
Fair — correctness is unaffected. I just meant the partial.size == 0 case is untested; a band with only fully-contained faces would cover it. Non-blocking.
Collaborator
Author
There was a problem hiding this comment.
I think it technically is tested for bands=[-90, 90] in test_conservative_full_sphere_conservation, but maybe we can add a test on a less trivial band.
On the other hand, this opened a rabbit hole that led to issue #1644, so we might want to fix that first.
erogluorhan
reviewed
Aug 7, 2026
erogluorhan
left a comment
erogluorhan
left a comment
Member
There was a problem hiding this comment.
Overall this looks good to me. Just a couple questions:
- Should we worry about the couple benchmarks that got worse? I believe, no, they shouldn't be directly related to the changes here, but am curious about your thoughts.
- Could you add test cases to cover the two latent bugs that this PR fixes?
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
4 participants
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.
This PR contains two post-accusphere optimizations, eliminating low-level hard materializations by using a vector-based masking strategy rather than individual conditionals, and reduced
zonal_meanpeakmem by building only the candidate faces instead of the whole grid.Partly addresses #1587
Closes #1650
Overview
PR Checklist
General
Testing
Documentation