Skip to content

add numba_math.py; optimize coordinates.py and _small_angle_of_2_vectors - #1674

Merged
Sevans711 merged 15 commits into
mainfrom
sevans/add-numba_math.py--optimize-coordinates.py
Aug 18, 2026
Merged

add numba_math.py; optimize coordinates.py and _small_angle_of_2_vectors#1674
Sevans711 merged 15 commits into
mainfrom
sevans/add-numba_math.py--optimize-coordinates.py

Conversation

@Sevans711

@Sevans711Sevans711 commented Aug 13, 2026

Copy link
Copy Markdown
Collaborator

Closes#1680 and closes#1681 (both are sub-issues of #1648).

Overview

(For #1680) adds uxarray/utils/numba_math.py to provide numba routines primarily for 3-vectors. Applying these directly to tuples to avoid constructing lots of tiny numpy arrays inside numba routines can cause huge speedup in some cases.

(For #1681) Optimizes numba routines in uxarray/grid/coordinates.py, and the _small_angle_of_2_vectors function (which is called directly by some routines in coordinates.py). As per discussion in #1648, this PR does not yet attempt to make these optimizations more broadly, but instead just focuses on improving a single area of the code.

After these changes, locally I am seeing grid.construct_face_centers('welzl') time decrease from ~14 s to ~150 ms, for a grid with 157k faces. ASV benchmarks similarly show huge improvement (>=20x speedup) there.

Sidenote: attempted to make a cleverer / more convenient interface, such as def _numba_add3(*vectors) which could add multiple vectors together, instead of just def _numba_add3(u, v) but that leads to extremely confusing numba errors in case of any type mismatch between values in tuples. Seems like a better idea to keep this simple (yet maybe slightly less convenient) interface to avoid those confusing errors.

Tiny expansion of PR scope: also fixes type-hinting for _lonlat_rad_to_xyz and _xyz_to_lonlat_rad.

PR Checklist

General

  • An issue is created and linked
  • Added appropriate labels (if your uxarray repo permissions allow it)
  • Filled out Overview and Expected Usage (if applicable) sections

Testing & Benchmarking

  • Adequate tests are created if there is new functionality
  • Tests are not too basic (such as simply calling a function and nothing else)
  • Tests cover all major paths in your new functions
  • If this PR could affect performance, ran ASV benchmarks and confirmed they show expected behavior (add a new benchmark if necessary)

Documentation

  • Docstrings have been added to all new functions
  • Docstrings have been updated with any function changes
  • [N/A] User (public) functions have been added to docs/api.rst
  • Internal (private) function names start with an underscore (_)

AI Disclosure

AI Usage: GitHub Copilot's inline code suggestions; discussions with Claude Sonnet and ChatGPT to help understand and debug numba routines and confusing numba errors.

  • I take responsibility for all AI-generated content in my PR.
  • I have tested all AI-generated content in my PR.

E.g. locally seeing grid.construct_face_centers('welzl') time decrease from 14 s to ~150 ms, for a grid with 157k faces.
@Sevans711Sevans711 added scalability Related to scalability & performance efforts run-benchmark Run ASV benchmark workflow labels Aug 13, 2026
@github-actions

github-actionsBot commented Aug 13, 2026

Copy link
Copy Markdown

ASV Benchmarking

Benchmark Comparison Results

Benchmarks that have improved:

ChangeBefore [ab39109]After [07a04fc]RatioBenchmark (Parameter)
-1.49±0.03s71.5±2ms0.05mpas_ocean.ConstructFaceLatLon.time_welzl('120km')
-92.7±0.6ms6.51±0.1ms0.07mpas_ocean.ConstructFaceLatLon.time_welzl('480km')

Benchmarks that have stayed the same:

ChangeBefore [ab39109]After [07a04fc]RatioBenchmark (Parameter)
122±10ms113±1ms0.92bench_connectivity.Connectivity.time_edge_face('120km')
7.13±0.09ms6.79±0.04ms0.95bench_connectivity.Connectivity.time_edge_face('480km')
115±3ms114±2ms1.00bench_connectivity.Connectivity.time_edge_node('120km')
6.59±0.3ms6.33±0.09ms0.96bench_connectivity.Connectivity.time_edge_node('480km')
116±1ms113±1ms0.98bench_connectivity.Connectivity.time_face_edge('120km')
6.66±0.2ms6.40±0.03ms0.96bench_connectivity.Connectivity.time_face_edge('480km')
468±20ms455±4ms0.97bench_connectivity.Connectivity.time_face_face('120km')
28.0±1ms27.2±0.2ms0.97bench_connectivity.Connectivity.time_face_face('480km')
48.8±7μs43.4±2μs~0.89bench_connectivity.Connectivity.time_face_node('120km')
49.5±4μs40.3±3μs~0.81bench_connectivity.Connectivity.time_face_node('480km')
286±40μs258±5μs~0.90bench_connectivity.Connectivity.time_n_nodes_per_face('120km')
211±20μs210±9μs0.99bench_connectivity.Connectivity.time_n_nodes_per_face('480km')
117±4ms116±0.5ms0.99bench_connectivity.Connectivity.time_node_edge('120km')
7.08±0.5ms6.77±0.3ms0.96bench_connectivity.Connectivity.time_node_edge('480km')
46.3±5ms51.0±0.9ms~1.10bench_connectivity.Connectivity.time_node_face('120km')
2.81±0.2ms2.58±0.01ms0.92bench_connectivity.Connectivity.time_node_face('480km')
5.17±1ms4.81±0.1ms0.93face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/mpas/QU/oQU480.231010.nc'))
1.69±0.3ms1.55±0.05ms0.92face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/scrip/outCSne8/outCSne8.nc'))
6.51±7s6.73±8ms~0.00face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/geoflow-small/grid.nc'))
1.12±0.03ms1.16±0.02ms1.04face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/quad-hexagon/grid.nc'))
57.3k57.3k1.00face_bounds.FaceBounds.track_nbytes_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/mpas/QU/oQU480.231010.nc'))
12.3k12.3k1.00face_bounds.FaceBounds.track_nbytes_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/scrip/outCSne8/outCSne8.nc'))
123k123k1.00face_bounds.FaceBounds.track_nbytes_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/geoflow-small/grid.nc'))
1281281.00face_bounds.FaceBounds.track_nbytes_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/quad-hexagon/grid.nc'))
1.27M1.27M1.00face_bounds.FaceBounds.track_nbytes_grid_with_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/mpas/QU/oQU480.231010.nc'))
50.1k50.1k1.00face_bounds.FaceBounds.track_nbytes_grid_with_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/scrip/outCSne8/outCSne8.nc'))
1.48M1.48M1.00face_bounds.FaceBounds.track_nbytes_grid_with_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/geoflow-small/grid.nc'))
7127121.00face_bounds.FaceBounds.track_nbytes_grid_with_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/quad-hexagon/grid.nc'))
2.04M2.03M1.00face_bounds.FaceBounds.track_peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/mpas/QU/oQU480.231010.nc'))
2.04M2.03M1.00face_bounds.FaceBounds.track_peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/scrip/outCSne8/outCSne8.nc'))
2.21M2.2M1.00face_bounds.FaceBounds.track_peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/geoflow-small/grid.nc'))
38.3k38.3k1.00face_bounds.FaceBounds.track_peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/quad-hexagon/grid.nc'))
337M336M1.00face_bounds.FaceBoundsColdStartRss.peakmem_open_and_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/mpas/QU/oQU480.231010.nc'))
367M368M1.00face_bounds.FaceBoundsColdStartRss.peakmem_open_and_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/scrip/outCSne8/outCSne8.nc'))
338M338M1.00face_bounds.FaceBoundsColdStartRss.peakmem_open_and_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/geoflow-small/grid.nc'))
337M337M1.00face_bounds.FaceBoundsColdStartRss.peakmem_open_and_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/quad-hexagon/grid.nc'))
671±40ns676±20ns1.01geometry_kernels.AccucrossKernels.time_accucross
1.38±0.02μs1.43±0.01μs1.03geometry_kernels.AccucrossKernels.time_accucross_pair
256±10ns245±4ns0.96geometry_kernels.EFTPrimitives.time_acc_sqrt_re
240±3ns250±8ns1.04geometry_kernels.EFTPrimitives.time_diff_of_products
216±5ns225±5ns1.04geometry_kernels.EFTPrimitives.time_two_prod
226±5ns221±10ns0.98geometry_kernels.EFTPrimitives.time_two_sum
792±9ns851±5ns1.07geometry_kernels.GCAConstLatIntersection.time_accux_constlat_kernel
576±30ns611±20ns1.06geometry_kernels.GCAConstLatIntersection.time_gca_const_lat_intersection
1.00±0.02μs1.05±0.02μs1.05geometry_kernels.GCAConstLatIntersection.time_try_gca_const_lat_intersection
866±30ns876±10ns1.01geometry_kernels.GCAGCAIntersection.time_accux_gca_kernel
796±20ns806±10ns1.01geometry_kernels.GCAGCAIntersection.time_gca_gca_intersection
1.16±0.02μs1.14±0.02μs0.98geometry_kernels.GCAGCAIntersection.time_try_gca_gca_intersection
24.0±0.7μs22.8±1μs0.95geometry_kernels.OrientPredicates.time_on_minor_arc
461±9ns436±10ns0.95geometry_kernels.OrientPredicates.time_orient3d_on_sphere
1.61±0.04ms1.57±0.04ms0.98geometry_samebody.SameBodyConstLat.time_accux_dispatch
681±10μs681±5μs1.00geometry_samebody.SameBodyConstLat.time_accux_kernel
1.11±0.01ms1.11±0.03ms1.00geometry_samebody.SameBodyConstLat.time_fp64_dispatch
90.2±3μs86.9±2μs0.96geometry_samebody.SameBodyConstLat.time_fp64_kernel
25.5±0.3ms25.4±0.7ms1.00geometry_samebody_gcagca.SameBodyGcaGca.time_accux_dispatch
5.91±0.1ms5.76±0.1ms0.97geometry_samebody_gcagca.SameBodyGcaGca.time_accux_kernel
21.9±0.1ms21.3±0.4ms0.97geometry_samebody_gcagca.SameBodyGcaGca.time_fp64_dispatch
3.59±0.1ms3.52±0.03ms0.98geometry_samebody_gcagca.SameBodyGcaGca.time_fp64_kernel
555±20ms532±20ms0.96import.Imports.timeraw_import_uxarray
293M293M1.00import.Imports.track_peakmem_import_uxarray
1.45±0.03ms1.48±0.02ms1.02mpas_ocean.CheckNorm.time_check_norm('120km')
1.15±0.02ms1.12±0.02ms0.98mpas_ocean.CheckNorm.time_check_norm('480km')
419±20ms406±10ms0.97mpas_ocean.ConnectivityConstruction.time_face_face_connectivity('120km')
25.1±1ms25.2±0.7ms1.01mpas_ocean.ConnectivityConstruction.time_face_face_connectivity('480km')
400±9μs402±10μs1.00mpas_ocean.ConnectivityConstruction.time_n_nodes_per_face('120km')
346±10μs330±6μs0.95mpas_ocean.ConnectivityConstruction.time_n_nodes_per_face('480km')
2.62±0.01ms2.68±0.08ms1.02mpas_ocean.ConstructFaceLatLon.time_cartesian_averaging('120km')
2.19±0.03ms2.13±0.07ms0.98mpas_ocean.ConstructFaceLatLon.time_cartesian_averaging('480km')
13.2±0.2ms12.9±0.3ms0.97mpas_ocean.ConstructTreeStructures.time_ball_tree('120km')
807±10μs799±40μs0.99mpas_ocean.ConstructTreeStructures.time_ball_tree('480km')
6.84±0.06ms6.64±0.06ms0.97mpas_ocean.ConstructTreeStructures.time_kd_tree('120km')
559±6μs535±20μs0.96mpas_ocean.ConstructTreeStructures.time_kd_tree('480km')
322±10ms334±9ms1.04mpas_ocean.CrossSections.time_const_lat('120km', 1)
166±3ms173±6ms1.04mpas_ocean.CrossSections.time_const_lat('120km', 2)
87.8±1ms86.7±5ms0.99mpas_ocean.CrossSections.time_const_lat('120km', 4)
276±7ms276±7ms1.00mpas_ocean.CrossSections.time_const_lat('480km', 1)
139±2ms144±2ms1.03mpas_ocean.CrossSections.time_const_lat('480km', 2)
71.5±2ms74.2±0.5ms1.04mpas_ocean.CrossSections.time_const_lat('480km', 4)
357M357M1.00mpas_ocean.CrossSectionsPeakMem.peakmem_const_lat('120km', 1)
356M357M1.00mpas_ocean.CrossSectionsPeakMem.peakmem_const_lat('120km', 2)
357M356M1.00mpas_ocean.CrossSectionsPeakMem.peakmem_const_lat('120km', 4)
340M342M1.01mpas_ocean.CrossSectionsPeakMem.peakmem_const_lat('480km', 1)
340M339M1.00mpas_ocean.CrossSectionsPeakMem.peakmem_const_lat('480km', 2)
340M340M1.00mpas_ocean.CrossSectionsPeakMem.peakmem_const_lat('480km', 4)
16.0±0.3ms16.6±0.3ms1.04mpas_ocean.DualMesh.time_dual_mesh_construction('120km')
1.78±0.07ms1.89±0.08ms1.06mpas_ocean.DualMesh.time_dual_mesh_construction('480km')
49.8±5ms48.5±0.4ms0.97mpas_ocean.FaceAreas.time_face_areas('120km')
5.23±0.4ms4.94±0.2ms0.94mpas_ocean.FaceAreas.time_face_areas('480km')
229k229k1.00mpas_ocean.FaceAreas.track_nbytes_face_areas('120km')
14.3k14.3k1.00mpas_ocean.FaceAreas.track_nbytes_face_areas('480km')
2.12M2.12M1.00mpas_ocean.FaceAreas.track_peakmem_face_areas('120km')
836k836k1.00mpas_ocean.FaceAreas.track_peakmem_face_areas('480km')
470±10ms480±4ms1.02mpas_ocean.GeoDataFrame.time_to_geodataframe('120km', False)
29.2±1ms31.8±4ms1.09mpas_ocean.GeoDataFrame.time_to_geodataframe('120km', True)
43.0±2ms41.5±0.8ms0.96mpas_ocean.GeoDataFrame.time_to_geodataframe('480km', False)
4.09±0.1ms3.89±0.05ms0.95mpas_ocean.GeoDataFrame.time_to_geodataframe('480km', True)
139±3ms135±4ms0.97mpas_ocean.Gradient.time_gradient('120km')
9.14±0.4ms9.37±0.5ms1.03mpas_ocean.Gradient.time_gradient('480km')
457k457k1.00mpas_ocean.Gradient.track_nbytes_gradient('120km')
28.7k28.7k1.00mpas_ocean.Gradient.track_nbytes_gradient('480km')
5.08M5.08M1.00mpas_ocean.Gradient.track_peakmem_gradient('120km')
328k328k1.00mpas_ocean.Gradient.track_peakmem_gradient('480km')
352M352M1.00mpas_ocean.GradientColdStartRss.peakmem_gradient('120km')
331M331M1.00mpas_ocean.GradientColdStartRss.peakmem_gradient('480km')
408±20μs487±10μs~1.20mpas_ocean.Integrate.time_integrate('120km')
18.4M18.4M1.00mpas_ocean.Integrate.track_nbytes_integrate('120km')
1.2M1.2M1.00mpas_ocean.Integrate.track_nbytes_integrate('480km')
99.9±0.9ms107±2ms1.07mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('120km', 'exclude')
101±1ms105±2ms1.04mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('120km', 'include')
102±0.6ms109±3ms1.07mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('120km', 'split')
8.18±0.2ms8.75±0.2ms1.07mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('480km', 'split')
241±8μs255±20μs1.06mpas_ocean.PointInPolygon.time_face_search_lonlat('120km')
203±1μs222±8μs1.09mpas_ocean.PointInPolygon.time_face_search_lonlat('480km')
206±5μs206±6μs1.00mpas_ocean.PointInPolygon.time_face_search_xyz('480km')
122±1ms126±0.6ms1.03mpas_ocean.RemapDownsample.time_bilinear_remapping
121±2ms124±0.5ms1.02mpas_ocean.RemapDownsample.time_inverse_distance_weighted_remapping
9.56±0.2ms10.3±0.3ms1.08mpas_ocean.RemapDownsample.time_nearest_neighbor_remapping
620±8ms615±8ms0.99mpas_ocean.RemapUpsample.time_bilinear_remapping
23.7±0.4ms23.8±0.4ms1.00mpas_ocean.RemapUpsample.time_inverse_distance_weighted_remapping
7.30±0.3ms7.73±0.5ms1.06mpas_ocean.RemapUpsample.time_nearest_neighbor_remapping
5.17±0.4ms5.66±0.09ms1.09mpas_ocean.ZonalAverage.time_zonal_average('120km')
2.61±0.2ms2.75±0.1ms1.05mpas_ocean.ZonalAverage.time_zonal_average('480km')
358M358M1.00mpas_ocean.ZonalAveragePeakMem.peakmem_zonal_average('120km')
341M341M1.00mpas_ocean.ZonalAveragePeakMem.peakmem_zonal_average('480km')
3.68±0.04ms3.99±0.3ms1.08quad_hexagon.QuadHexagon.time_open_dataset
3.13±0.02ms3.14±0.2ms1.00quad_hexagon.QuadHexagon.time_open_grid
4084081.00quad_hexagon.QuadHexagon.track_nbytes_open_dataset
3923921.00quad_hexagon.QuadHexagon.track_nbytes_open_grid
73.7k73.5k1.00quad_hexagon.QuadHexagon.track_peakmem_open_dataset
73k72.8k1.00quad_hexagon.QuadHexagon.track_peakmem_open_grid

Benchmarks that have got worse:

ChangeBefore [ab39109]After [07a04fc]RatioBenchmark (Parameter)
+258±9μs285±7μs1.1mpas_ocean.HoleEdgeIndices.time_construct_hole_edge_indices('120km')
+131±8μs151±8μs1.15mpas_ocean.HoleEdgeIndices.time_construct_hole_edge_indices('480km')
+307±7μs348±20μs1.13mpas_ocean.Integrate.time_integrate('480km')
+7.90±0.05ms9.08±0.4ms1.15mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('480km', 'exclude')
+7.85±0.2ms8.92±0.1ms1.14mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('480km', 'include')
+209±7μs318±50μs1.52mpas_ocean.PointInPolygon.time_face_search_xyz('120km')

@Sevans711
Sevans711 marked this pull request as ready for review August 14, 2026 15:25
@Sevans711Sevans711 added run-benchmark Run ASV benchmark workflow and removed run-benchmark Run ASV benchmark workflow labels Aug 17, 2026
Comment threaduxarray/grid/coordinates.py Outdated
Comment threaduxarray/grid/coordinates.py Outdated
@cmdupuis3cmdupuis3 added run-benchmark Run ASV benchmark workflow and removed run-benchmark Run ASV benchmark workflow labels Aug 17, 2026
@cmdupuis3

cmdupuis3 commented Aug 17, 2026

Copy link
Copy Markdown
Collaborator

Also, I think if we're going to go this route with numba_math, we should probably have a standard vector type for the whole repo, and then for these binops we can use operator overloads like @overload(operator.sub) so we don't have to explicitly call a long-named function when we just want to subtract things. Either way seems like a big change to the repo philosophy though.

@cmdupuis3

Copy link
Copy Markdown
Collaborator

Ope, wait, there's an issue with operator overloads. If you disable the numba jit, the decorator skips and Python goes back to concatenating with + instead of adding.

The joys of Python...

@Sevans711

Copy link
Copy Markdown
CollaboratorAuthor

Thank you @cmdupuis3 for reviewing!

Also, I think if we're going to go this route with numba_math, we should probably have a standard vector type for the whole repo, and then for these binops we can use operator overloads like @overload(operator.sub) so we don't have to explicitly call a long-named function when we just want to subtract things. Either way seems like a big change to the repo philosophy though.

Ope, wait, there's an issue with operator overloads. If you disable the numba jit, the decorator skips and Python goes back to concatenating with + instead of adding.

The joys of Python...

Yeahhh, while Python objects would be happy to support a convenient interface with overloads like that, it seems to me like for performance-optimized numba it might be better to just keep the clunkier interface. Numba compilations make things harder to reason about, so keeping it simple seems to be the way to go here.

(E.g. I also originally tried def _numba_add3(*vectors, scalar=None) as a way to allow adding any number of 3-vectors plus maybe a scalar if you want, but this led to extremely confusing behavior, such as _numba_add3((1,2,3),scalar=10) crashing while _numba_add3((1,2,3),10) unexpectedly ran without complaint and returned (11,12,13); and _numba_add3((1,2,3),(4,5,6),(7,8,9)) working okay but _numba_add3((1,2,3),(4,5,6),(7,8,9.0)) crashing with an extremely confusing numba typing error)

To me it doesn't feel like a change to the whole repo philosophy, only a change to the philosophy of numba-compiled functions (wherein this directly causes huge speedup). Outside of njit-decorated functions, I believe uxarray should certainly still use numpy/dask/xarray arrays for dealing with vectors!

@erogluorhanerogluorhan left a comment

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

This looks good to me, so I approve it; just a single comment below though.

Comment threaduxarray/grid/coordinates.py Outdated
@Sevans711
Sevans711 merged commit b674956 into mainAug 18, 2026
13 checks passed
@Sevans711
Sevans711 deleted the sevans/add-numba_math.py--optimize-coordinates.py branch August 18, 2026 20:38
Sign up for freeto join this conversation on GitHub. Already have an account? Sign in to comment

Labels

run-benchmarkRun ASV benchmark workflowscalabilityRelated to scalability & performance efforts

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Avoid tiny numpy arrays in numba: optimize coordinates.py Avoid tiny numpy arrays in numba: create numba helper routines file

3 participants

@Sevans711@cmdupuis3@erogluorhan