Conversation
Match nodes in Cartesian space rather than the lon/lat plane so pole and antimeridian nodes are recognized as the same point, and store the resulting polar faces as triangles instead of quads with a repeated corner.
Canonicalize duplicate node indices in the face-node connectivity before building the dual, so grids with repeated nodes produce a correct dual instead of being rejected. Also vectorize the duplicate lookup and deprecate the now redundant check_duplicate_nodes argument.
Sevans711
left a comment
Sevans711
left a comment
There was a problem hiding this comment.
Mostly looks like a good fix! I think it still needs a little bit of extra work, and I left some inline comments accordingly.
| """Map duplicate node indices to the first index with the same coordinates.""" | ||
| node_coordinates = np.column_stack((grid.node_lon.values, grid.node_lat.values)) | ||
| _, first_indices, inverse_indices = np.unique( | ||
| node_coordinates, axis=0, return_index=True, return_inverse=True |
There was a problem hiding this comment.
Is "exact equality" the correct way to go here? My intuition originally was that there should probably be some sort of tolerance here, e.g. if values agree to within 1e-12, they are probably the same node, right?
There was a problem hiding this comment.
In fact, looking back at the original issue thread, it looks like you were the person who originally suggested there be a tolerance in the first place! So, I would actually now assert that exact inequality is not the desired implementation, there should be a tolerance, as clarified in thread of #865.
There was a problem hiding this comment.
Fixed — matching goes through _coincident_node_canonical_indices, which is a KDTree query within ERROR_TOLERANCE on the sphere, not exact equality.
There was a problem hiding this comment.
This fix looks reasonable. One more nitpick here though: your comment on the original thread #865 (comment) suggested you wanted the user to be able to specify the tolerance.
Right now, there is no way to specify the "tolerance for merging duplicate nodes" from any public-facing functions. To specify it as a larger value (e.g. 1e-4) a user would need to call _merge_coincident_grid_ds_nodes directly on grid_ds before passing to Grid.__init__, while there is presumably no way to properly specify it as a smaller value (e.g. 1e-12) because the _merge_coincident_grid_ds_nodes call inside Grid.__init__ would always rerun the process with a tolerance of ERROR_TOLERANCE=1e-8.
Should this option be added somewhere (e.g. something like a duplicates_tolerance parameter in Grid.__init__) or is "allow the user to specify the tolerance" no longer desirable?
There was a problem hiding this comment.
Can be a good edition, but not in this PR. I'll open a separate issue for a user-facing tolerance.
Adding this means threading it through open_grid and the from_* constructors..
There was a problem hiding this comment.
(Keeping this unresolved as a reminder to open that issue. Feel free to mark as resolved as soon as you opened that issue!)
| np.arange(grid.n_node, dtype=INT_DTYPE) != first_indices[inverse_indices] | ||
| ) | ||
| return { | ||
| INT_DTYPE(index): INT_DTYPE(first_indices[inverse_indices[index]]) |
There was a problem hiding this comment.
I would guess that creating a dict here is extremely inefficient… maybe that doesn't matter if there are only ever a tiny number of duplicate nodes. Not necessarily blocking, but have you looked into how long this takes to run for any larger grids containing duplicate nodes? (Or, do you expect a very limited number of duplicate nodes in most cases? I would guess it is probably not worth worrying about if there are less than ~1000 duplicates or so.)
There was a problem hiding this comment.
The dict only holds duplicates, not every node, so it stays small. The real cost was _check_duplicate_nodes_indices looping over every face in Python; that is now one np.isin over the connectivity array.
There was a problem hiding this comment.
I think that makes sense. Added the run-benchmarks label because this PR adds a call to _merge_coincident_grid_ds_nodes into Grid.__init__ so it might affect uxarray performance. I believe my original comment here is resolved, but I am keeping this thread open for now as a reminder to myself to consider the benchmarking results before approving this PR.
There was a problem hiding this comment.
The cache wasn't the whole story - dual mesh is still ~1.6x slower because get_dual builds a new Grid, which re-runs the coincident merge on the dual. Fixing that next, so don't go by the current benchmark numbers.
There was a problem hiding this comment.
That fix definitely helped, but it still looks like >40% slowdown in dual_mesh construction. Any insights?
Though, now that the performance improvements are expanding significantly beyond the original lines I flagged here, I think it might make more sense to continue the benchmarking discussion in the main thread below instead. Feel free to mark this as resolved once you read it, I will leave a comment soon during re-review.
|
Actually, apologies for not including this during the original review, comment but one more thought: does this actually fully close the original issue? The issue writeup makes it sound to me like duplicate nodes should be handled immediately upon constructing the grid, not just during one functionality (get_dual). Is there a reason that duplicate nodes should be handled only during get_dual, instead of immediately? (Do all other current/planned functions work properly regardless of whether there are duplicate nodes?) |
float32 input (e.g. real climate datasets) silently ran the whole xyz/tolerance pipeline at float32 precision, causing pole/antimeridian merges to fail or merge only partially.
Match nodes in Cartesian space rather than the lon/lat plane so pole and antimeridian nodes are recognized as the same point, and store the resulting polar faces as triangles instead of quads with a repeated corner.
float32 input (e.g. real climate datasets) silently ran the whole xyz/tolerance pipeline at float32 precision, causing pole/antimeridian merges to fail or merge only partially.
Extend #865's fix beyond the dual mesh: canonicalize duplicate/coincident node indices in connectivity for every Grid construction path, not just construct_dual. Detection is now tolerance-based (unit-sphere chordal distance) instead of exact lon/lat match, so pole-degenerate duplicates are also caught. Node coordinate arrays are left untouched by design; only connectivity is remapped to canonical indices, with any resulting repeated face corners collapsed.
construct_dual no longer needs its own per-call duplicate detection and remap, and get_dual() no longer needs to hard-gate on duplicate node indices, since Grid construction now canonicalizes them structurally before any of this code runs.
Since duplicate node coordinates are intentionally left unreferenced by connectivity, a node KDTree/BallTree built over the raw coordinate array could select an index no face actually points to, silently returning empty or wrong nearest-neighbor results. Build the "nodes" tree only over live (referenced) indices and translate query results back to original index space.
polars' unique() with maintain_order unset does not guarantee row order across runs, so the node index assigned to a given corner coordinate could vary between reads of the same file. This is normally harmless, but it made canonical-node selection for coincident duplicates (e.g. pole points with differing longitude) flaky from run to run.
test_dual_duplicate: validate() now succeeds since connectivity is fully canonicalized (duplicate coordinates remain by design, but nothing references a dead index anymore). test_grid_nn_subset: max valid k for a node search is now bounded by the live node count, not raw node count. test_to_geodataframe_preserves_antimeridian_faces: pole-coincident corners with differing longitude are now also merged, shifting the antimeridian face count.
Merging pole-adjacent duplicate nodes was collapsing each face's own locally-meaningful longitude at the pole into one arbitrary canonical value, which corrupted lat/lon bounds and broke zonal weight computation for cube-sphere grids near the poles.
Sevans711
left a comment
Sevans711
left a comment
There was a problem hiding this comment.
Suggestion: please change the title of this PR to reflect the full scope. Duplicate nodes are now handled directly whenever constructing a Grid, not just when constructing the dual mesh. I was confused about why _check_duplicate_nodes_indices checks were removed despite construct_dual() being unchanged. I think the reason is because all Grid objects are now guaranteed to not have duplicate nodes.
I tried leaving a full review but I kept getting the feeling that something weird was happening, a nagging feeling like "hey I think I've looked at this code before and left an inline comment, why am I reviewing it again?" Then I realized many of the changes here are also in #1690 which hasn't been merged yet. That led to duplicating review work and will probably lead to needing to apply fixes multiple times. I'm not sure the cleanest way forwards at this point, but I might suggest waiting to merge this until after #1690 gets merged. At least, I will want to do another close review of this PR after that PR merges, because I lost track of which things I already looked at closely and which I need to consider again.
…des' into rajeeja/coincident-nodes # Conflicts: # test/io/test_structured.py # uxarray/core/dataarray.py # uxarray/core/dataset.py # uxarray/grid/grid.py # uxarray/io/_structured.py
"Coincident nodes" is the domain term used by TempestRemap, MOAB, and our own validation module; "dedupe" reads like a dataframe operation. Also note in the docstring that we keep the redundant coordinates and only remap connectivity, where TempestRemap and MOAB delete and renumber.
_coincident_node_canonical_indices exempts pole nodes from merging so each face touching a pole keeps its own longitude. The mask used np.isclose(|z|, 1.0, atol=tolerance), which is wrong twice: tolerance is a chord radius everywhere else in the function, and np.isclose's default rtol=1e-5 swamped atol entirely. The carve-out spanned 1 - |z| <= 1.001e-5, a chord of 4.5e-3 or ~28 km on Earth, so any node within 28 km of a pole was silently exempt from merging. For a point at colatitude t, 1 - |z| = 1 - cos(t) = chord**2 / 2, so the correct cap is tolerance**2 / 2 with rtol pinned to zero.
construct_dual reads node_face_connectivity with no duplicate handling, so a face still referencing a coincident duplicate index yields a degenerate dual face instead of an error. Merging at construction makes that unreachable today, so the guard is a no-op safety net; check_duplicate_nodes is now deprecated and ignored, since the check always runs.
The pole/seam merge in _read_structured_grid changes the ersstv5 grid (16290 nodes to 16200) and turns pole quads into triangles, so the near-pole plots were stale. Outputs re-executed; the stray OMP stderr line is dropped.
|
Check out this pull request on See visual diffs & provide feedback on Jupyter Notebooks. Powered by ReviewNB |
The notebook's cell sources are identical to main; only the outputs differed, and those differences were an xarray coordinate-ordering repr change and a stray OMP info line. Keeping the re-run added ~1700 lines of diff with no content change.
Coincident-node merging now runs on every Grid construction, so the cost is paid by grids that have no duplicates at all. Two points within a chord of ERROR_TOLERANCE differ by at most that in x, so in x-sorted order every consecutive gap between them is also within tolerance; a point whose sorted neighbours are both further away cannot be coincident with anything. Filtering on that is exact and drops the tree build entirely for a clean grid: 1.23s -> 0.28s over 2M nodes.
rajeeja
changed the title
ASV BenchmarkingBenchmark Comparison ResultsBenchmarks that have stayed the same:
Benchmarks that have got worse:
|
Sevans711
left a comment
Sevans711
left a comment
There was a problem hiding this comment.
Making good progress! Still a few inline comments to resolve, see below and also I left some follow-up comments on previous inline comment threads.
Also, looking at ASV benchmark results, it looks like quite a few benchmarks have slowed down significantly, especially dual mesh construction being ~50% slower. Why would dual mesh construction slow down so much, if the search for duplicates primarily occurs when building the original grid? Tagging @cmdupuis3 to maybe comment on the benchmarks if you have additional insights. See also #1710 as a sanity check that the benchmark changes represent real changes to performance, rather than quirks of the benchmarking suite.
Thanks, the dual mesh slowdown was the duplicate check running on every get_dual call and rescanning every node each time. It's cached on the grid now. |
get_dual constructs a new Grid from this grid's face centers, which re-ran the coincident node search on every call and made dual mesh construction ~1.6x slower. The dual's connectivity is already canonical, so the merge is skipped there via a new internal Grid flag.
…ason through corner collapsing
Sevans711
left a comment
Sevans711
left a comment
There was a problem hiding this comment.
One minor inline comment. Aside from that, my only remaining concern is about performance.
Thanks, the dual mesh slowdown was the duplicate check running on every get_dual call and rescanning every node each time. It's cached on the grid now.
This seems to have helped, but the latest benchmarks still shows a >40% slowdown. Additionally, the slowdown is worse for the larger grid (120km resolution reported 45% slowdown, 480km resolution reported 34% slowdown), which suggests the slowdown might be much worse for larger, production-scale grids.
This led me to trying a larger grid as a sanity check. I have a 15km MPAS grid on my local machine (it's the dyamond 15km grid mentioned in the benchmarks suite) so I used that. Running ux.open_grid() a few times, I'm seeing it takes roughly 600 ms on main, and roughly 3.6 s on this branch. That's a 6x slowdown! Considering the possibility for nonlinear scaling, and that people often use even better resolution than 15km (e.g. 7.5km or 3.75km), this may be a significant performance hit for workflows with large grids.
Any ideas for how to improve this? I suspect that merging a 6x slowdown into open_grid() might not be worth always enforcing this extra validation, as many grids do not have any duplicates at all. (If there's really no way to improve the performance, maybe this check during Grid.__init__ should be disabled by default, with a flag to enable it?)
cc: @cmdupuis3
| for index in np.flatnonzero(canonical != np.arange(n_node, dtype=INT_DTYPE)) | ||
| } | ||
| if not duplicate_node_map: | ||
| return grid_ds |
There was a problem hiding this comment.
Minor nit: could you indent the subsequent lines to make it clear that the logic here looks like this:
if no duplicates exist:
no need to handle any duplicates
else:
handle duplicates
Alternatively, could swap the order if you want:
if not duplicate_node_map:
handle duplicates
else:
return grid_ds
With the current code, I had to scroll around a bit to find this line and understand that is what is happening.
Alternatively, if you added a small comment in this function's docstring to emphasize "if there are no duplicates, returns grid_ds, unchanged", that would be sufficient. (Or, could add a small comment below, or above when defining _NODE_INDEX_CONNECTIVITY_TO_REMAP and _DERIVED_CONNECTIVITY_TO_INVALIDATE, to clarify that these things only get applied if actually detected any duplicates.)
There was a problem hiding this comment.
To add, is there a way to say that some file formats will be guaranteed to not have duplicate nodes? If we know when we can skip this completely, that would be a nice thing to have.
3 participants
Closes #865
geoflow-smallpreviously could not produce a dual at all and now yields 3803 faces; the test asserts this and fails on main._find_duplicate_nodesis vectorized withnp.uniqueinstead of a per-node dict.Grid.get_dual(check_duplicate_nodes=...)is now ignored and deprecated rather than removed, so existing callers keep working.UxDataArray.get_dualandUxDataset.get_dualstill raiseGridInvalidError, since node-centered data cannot be remapped onto a merged node set.(lon, lat), such as poles and the antimeridian; that isGrid.from_structureddoes not merge coincident pole and antimeridian nodes #1689 / Merge coincident pole and antimeridian nodes in structured grids #1690.