From a8019c0ed43da3a316f788ac91fe5afba3ef491d Mon Sep 17 00:00:00 2001 From: cmdupuis3 Date: Wed, 19 Aug 2026 11:59:45 -0500 Subject: [PATCH 1/6] Remove per-edge allocations from gca_gca_intersection Add scalar-argument L1/L2 siblings (_accux_gca_scalar, _try_gca_gca_intersection_scalar) and rewire gca_gca_intersection to them, cutting 4 heap allocations per edge to 1. Bit-identical output; ~2.7x on the function in isolation, ~1.24x on serial face bounds. --- uxarray/grid/intersections.py | 175 +++++++++++++++++++++++++++++----- 1 file changed, 152 insertions(+), 23 deletions(-) diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py index 1e6100854..4313f8d84 100644 --- a/uxarray/grid/intersections.py +++ b/uxarray/grid/intersections.py @@ -313,6 +313,55 @@ def _accux_gca(w0, w1, v0, v1): return pos, neg +@njit(cache=True, inline="always", error_model="numpy") +def _accux_gca_scalar( + w00, w01, w02, w10, w11, w12, v00, v01, v02, v10, v11, v12 +): + """Scalar-argument form of :func:`_accux_gca`: returns the six components. + + Takes the twelve endpoint components directly and returns the candidate + components as scalars, so hot loops pay no heap allocation per edge. The + array form allocates two ``(3,)`` arrays per call; in the face-bounds path + that is two of the four allocations ``gca_gca_intersection`` used to make + for every edge of every face. + + Arithmetic is identical to :func:`_accux_gca`, operation for operation, so + results are bit-for-bit the same. + """ + n1x_hi, n1y_hi, n1z_hi, n1x_lo, n1y_lo, n1z_lo = accucross( + w00, w01, w02, w10, w11, w12 + ) + n2x_hi, n2y_hi, n2z_hi, n2x_lo, n2y_lo, n2z_lo = accucross( + v00, v01, v02, v10, v11, v12 + ) + vx_hi, vy_hi, vz_hi, vx_lo, vy_lo, vz_lo = accucross_pair( + n1x_hi, + n1y_hi, + n1z_hi, + n1x_lo, + n1y_lo, + n1z_lo, + n2x_hi, + n2y_hi, + n2z_hi, + n2x_lo, + n2y_lo, + n2z_lo, + ) + vx = vx_hi + vx_lo + vy = vy_hi + vy_lo + vz = vz_hi + vz_lo + sum_hi, sum_lo = _sum_of_squares_c((vx_hi, vy_hi, vz_hi), (vx_lo, vy_lo, vz_lo)) + vn, _ = acc_sqrt_re(sum_hi, sum_lo) + # vn==0 (coplanar arcs) yields inf via IEEE division under error_model="numpy", + # so the candidates become non-finite and the status layer masks them out. + inv = 1.0 / vn + pos_x = vx * inv + pos_y = vy * inv + pos_z = vz * inv + return pos_x, pos_y, pos_z, -pos_x, -pos_y, -pos_z + + @njit(cache=True, error_model="numpy") def _try_gca_gca_intersection(w0, w1, v0, v1): """Select the valid great-circle intersection and report a status code. @@ -361,6 +410,55 @@ def _try_gca_gca_intersection(w0, w1, v0, v1): return point, status, pos, neg +@njit(cache=True, inline="always", error_model="numpy") +def _try_gca_gca_intersection_scalar( + w00, w01, w02, w10, w11, w12, v00, v01, v02, v10, v11, v12 +): + """Scalar-argument form of :func:`_try_gca_gca_intersection`. + + Same mask arithmetic and the same status codes, but allocation-free: no + ``(3,)`` candidate arrays and no ``point`` array. Returns the selected point + components, the status code, and both candidates as scalars. + + Note the selected point is still formed by multiply-add masking, exactly as + in the array form. That is safe *only because* the caller branches on + ``status``: when both candidates are non-finite the masks are zero, + ``0.0 * nan`` makes the selected point ``nan``, and ``status == 2`` routes + the caller away from it. Do not replace the caller's status branch with mask + arithmetic -- a zero mask propagates a non-finite discarded operand rather + than discarding it. + """ + px, py, pz, ngx, ngy, ngz = _accux_gca_scalar( + w00, w01, w02, w10, w11, w12, v00, v01, v02, v10, v11, v12 + ) + + pos_fin = ( + int(math.isfinite(px)) * int(math.isfinite(py)) * int(math.isfinite(pz)) + ) + neg_fin = ( + int(math.isfinite(ngx)) * int(math.isfinite(ngy)) * int(math.isfinite(ngz)) + ) + pos_on_a = pos_fin * _on_minor_arc_xyz(px, py, pz, w00, w01, w02, w10, w11, w12) + pos_on_b = pos_fin * _on_minor_arc_xyz(px, py, pz, v00, v01, v02, v10, v11, v12) + neg_on_a = neg_fin * _on_minor_arc_xyz(ngx, ngy, ngz, w00, w01, w02, w10, w11, w12) + neg_on_b = neg_fin * _on_minor_arc_xyz(ngx, ngy, ngz, v00, v01, v02, v10, v11, v12) + + pos_valid = pos_fin * pos_on_a * pos_on_b + neg_valid = neg_fin * neg_on_a * neg_on_b + + pos_mask = pos_valid * (1 - neg_valid) + neg_mask = neg_valid * (1 - pos_valid) + + point_x = pos_mask * px + neg_mask * ngx + point_y = pos_mask * py + neg_mask * ngy + point_z = pos_mask * pz + neg_mask * ngz + + both = pos_valid * neg_valid + none = (1 - pos_valid) * (1 - neg_valid) + status = both + none * 2 + return point_x, point_y, point_z, status, px, py, pz, ngx, ngy, ngz + + @njit(cache=True, error_model="numpy") def gca_gca_intersection(gca_a_xyz, gca_b_xyz): """Return the intersection points of two great-circle arcs. @@ -396,40 +494,71 @@ def gca_gca_intersection(gca_a_xyz, gca_b_xyz): if gca_a_xyz.shape[1] != 3 or gca_b_xyz.shape[1] != 3: raise DimensionError("The two GCAs must be in the cartesian [x, y, z] format") - w0 = gca_a_xyz[0] - w1 = gca_a_xyz[1] - v0 = gca_b_xyz[0] - v1 = gca_b_xyz[1] - - point, status, pos, neg = _try_gca_gca_intersection(w0, w1, v0, v1) + # Unpack to scalars and run the allocation-free scalar chain. The array + # forms (_accux_gca, _try_gca_gca_intersection) allocated four (3,)/(2,3) + # arrays per call -- pos, neg, point, res -- which dominated this function + # in the face-bounds path, where it runs once per edge of every face. Only + # the (2, 3) result array remains, because the public return type is an + # array. The arithmetic is unchanged, so results are bit-for-bit identical. + w00 = gca_a_xyz[0, 0] + w01 = gca_a_xyz[0, 1] + w02 = gca_a_xyz[0, 2] + w10 = gca_a_xyz[1, 0] + w11 = gca_a_xyz[1, 1] + w12 = gca_a_xyz[1, 2] + v00 = gca_b_xyz[0, 0] + v01 = gca_b_xyz[0, 1] + v02 = gca_b_xyz[0, 2] + v10 = gca_b_xyz[1, 0] + v11 = gca_b_xyz[1, 1] + v12 = gca_b_xyz[1, 2] + + ( + point_x, + point_y, + point_z, + status, + pos_x, + pos_y, + pos_z, + neg_x, + neg_y, + neg_z, + ) = _try_gca_gca_intersection_scalar( + w00, w01, w02, w10, w11, w12, v00, v01, v02, v10, v11, v12 + ) res = np.empty((2, 3)) count = 0 + # The branch on ``status`` is load-bearing, not stylistic: when no candidate + # is valid both are non-finite, and the masked ``point_*`` above is nan. A + # mask-arithmetic rewrite of this dispatch would propagate that nan into the + # coplanar result (0.0 * nan == nan) instead of discarding it. if status == 0: - res[0, 0] = point[0] - res[0, 1] = point[1] - res[0, 2] = point[2] + res[0, 0] = point_x + res[0, 1] = point_y + res[0, 2] = point_z count = 1 elif status == 1: - res[0, 0] = pos[0] - res[0, 1] = pos[1] - res[0, 2] = pos[2] - res[1, 0] = neg[0] - res[1, 1] = neg[1] - res[1, 2] = neg[2] + res[0, 0] = pos_x + res[0, 1] = pos_y + res[0, 2] = pos_z + res[1, 0] = neg_x + res[1, 1] = neg_y + res[1, 2] = neg_z count = 2 else: # status == 2: no candidate on both arcs. # Check for coplanar overlap (shared endpoints) outside the kernel. - if on_minor_arc(v0, w0, w1): - res[count, 0] = v0[0] - res[count, 1] = v0[1] - res[count, 2] = v0[2] + if _on_minor_arc_xyz(v00, v01, v02, w00, w01, w02, w10, w11, w12): + res[count, 0] = v00 + res[count, 1] = v01 + res[count, 2] = v02 count += 1 - if on_minor_arc(v1, w0, w1): - res[count, 0] = v1[0] - res[count, 1] = v1[1] - res[count, 2] = v1[2] + if _on_minor_arc_xyz(v10, v11, v12, w00, w01, w02, w10, w11, w12): + res[count, 0] = v10 + res[count, 1] = v11 + res[count, 2] = v12 count += 1 return res[:count] From 1d6c3ea42e413144f7cef1000f7c3c1fa7898793 Mon Sep 17 00:00:00 2001 From: cmdupuis3 Date: Wed, 19 Aug 2026 12:03:18 -0500 Subject: [PATCH 2/6] Point GCAGCAIntersection L1/L2 benchmarks at the scalar kernels gca_gca_intersection no longer calls the array-form _accux_gca / _try_gca_gca_intersection, so benchmarking them no longer reflects the dispatcher's actual cost. Point at the scalar siblings instead, skipping gracefully (via skip_benchmark_if) on commits that predate them. --- benchmarks/geometry_kernels.py | 51 ++++++++++++++++++++++++---------- 1 file changed, 37 insertions(+), 14 deletions(-) diff --git a/benchmarks/geometry_kernels.py b/benchmarks/geometry_kernels.py index 08aea236e..5750bed46 100644 --- a/benchmarks/geometry_kernels.py +++ b/benchmarks/geometry_kernels.py @@ -14,6 +14,7 @@ """ import numpy as np +from asv_runner.benchmarks.mark import skip_benchmark_if def _unit(v): @@ -142,33 +143,55 @@ def time_on_minor_arc(self): self.on_minor_arc(_V0, _W0, _W1) +# ``gca_gca_intersection`` (Layer 3) calls scalar-argument L1/L2 kernels +# (``_accux_gca_scalar``, ``_try_gca_gca_intersection_scalar``) rather than the +# array-returning ``_accux_gca``/``_try_gca_gca_intersection`` that predate them +# -- the array forms each allocated a ``(3,)`` (or ``(2,3)``) result per call, +# which no longer happens on the hot path. Benchmarking the array forms would +# time functions the dispatcher does not call, so they are skipped -- not +# measured -- on commits before the scalar kernels existed. +try: + from uxarray.grid.intersections import ( # noqa: F401 + _accux_gca_scalar, + _try_gca_gca_intersection_scalar, + ) + + _HAS_SCALAR_GCA_KERNELS = True +except ImportError: + _HAS_SCALAR_GCA_KERNELS = False + + class GCAGCAIntersection: """Benchmark all three layers of the GCA-GCA intersection stack.""" def setup(self): - from uxarray.grid.intersections import ( - _accux_gca, - _try_gca_gca_intersection, - gca_gca_intersection, - ) + from uxarray.grid.intersections import gca_gca_intersection - self._accux_gca = _accux_gca - self._try_gca_gca_intersection = _try_gca_gca_intersection self.gca_gca_intersection = gca_gca_intersection - self.gca_a = np.stack([_W0, _W1]) self.gca_b = np.stack([_V0, _V1]) - _accux_gca(_W0, _W1, _V0, _V1) - _try_gca_gca_intersection(_W0, _W1, _V0, _V1) gca_gca_intersection(self.gca_a, self.gca_b) + if _HAS_SCALAR_GCA_KERNELS: + from uxarray.grid.intersections import ( + _accux_gca_scalar, + _try_gca_gca_intersection_scalar, + ) + + self._accux_gca_scalar = _accux_gca_scalar + self._try_gca_gca_intersection_scalar = _try_gca_gca_intersection_scalar + _accux_gca_scalar(*_W0, *_W1, *_V0, *_V1) + _try_gca_gca_intersection_scalar(*_W0, *_W1, *_V0, *_V1) + + @skip_benchmark_if(not _HAS_SCALAR_GCA_KERNELS) def time_accux_gca_kernel(self): - """Layer 1: pure numerical kernel.""" - self._accux_gca(_W0, _W1, _V0, _V1) + """Layer 1: pure numerical kernel (scalar form; allocation-free).""" + self._accux_gca_scalar(*_W0, *_W1, *_V0, *_V1) + @skip_benchmark_if(not _HAS_SCALAR_GCA_KERNELS) def time_try_gca_gca_intersection(self): - """Layer 2: batch/status layer.""" - self._try_gca_gca_intersection(_W0, _W1, _V0, _V1) + """Layer 2: batch/status layer (scalar form; allocation-free).""" + self._try_gca_gca_intersection_scalar(*_W0, *_W1, *_V0, *_V1) def time_gca_gca_intersection(self): """Layer 3: dispatcher (full public API).""" From ee56b2c7df40d23e454c9a193041e674e3fabf3f Mon Sep 17 00:00:00 2001 From: cmdupuis3 Date: Wed, 19 Aug 2026 14:27:56 -0500 Subject: [PATCH 3/6] scalarized gca_gca comments --- uxarray/grid/intersections.py | 56 ++++++++++++++++++++++------------- 1 file changed, 35 insertions(+), 21 deletions(-) diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py index 4313f8d84..1a8818bd4 100644 --- a/uxarray/grid/intersections.py +++ b/uxarray/grid/intersections.py @@ -319,14 +319,24 @@ def _accux_gca_scalar( ): """Scalar-argument form of :func:`_accux_gca`: returns the six components. - Takes the twelve endpoint components directly and returns the candidate - components as scalars, so hot loops pay no heap allocation per edge. The - array form allocates two ``(3,)`` arrays per call; in the face-bounds path - that is two of the four allocations ``gca_gca_intersection`` used to make - for every edge of every face. - - Arithmetic is identical to :func:`_accux_gca`, operation for operation, so - results are bit-for-bit the same. + Compute the candidate intersection points of two great-circle arcs. + + Pure numerical kernel (mirrors AccuSphGeom ``accux_gca``). + + Computes the two antipodal candidate intersection points of the great-circle + arcs w0-w1 and v0-v1. No branching, no validity filtering. + + Parameters + ---------- + w00, w01, w02, w10, w11, w12 : float + Cartesian endpoints of the first arc. + v00, v01, v02, v10, v11, v12 : float + Cartesian endpoints of the second arc. + + Returns + ------- + pos_x, pos_y, pos_z, neg_x, neg_y, neg_z : float + Two antipodal candidate unit vectors. """ n1x_hi, n1y_hi, n1z_hi, n1x_lo, n1y_lo, n1z_lo = accucross( w00, w01, w02, w10, w11, w12 @@ -416,27 +426,31 @@ def _try_gca_gca_intersection_scalar( ): """Scalar-argument form of :func:`_try_gca_gca_intersection`. - Same mask arithmetic and the same status codes, but allocation-free: no - ``(3,)`` candidate arrays and no ``point`` array. Returns the selected point - components, the status code, and both candidates as scalars. - - Note the selected point is still formed by multiply-add masking, exactly as - in the array form. That is safe *only because* the caller branches on - ``status``: when both candidates are non-finite the masks are zero, - ``0.0 * nan`` makes the selected point ``nan``, and ``status == 2`` routes - the caller away from it. Do not replace the caller's status branch with mask - arithmetic -- a zero mask propagates a non-finite discarded operand rather - than discarding it. + Select the valid great-circle intersection and report a status code. + + Batch/status layer (mirrors AccuSphGeom ``try_gca_gca_intersection``). + + Calls the pure numerical kernel, applies integer mask arithmetic to determine + validity, selects the output point without if/else branching in the hot path. + + Status codes mirror AccuSphGeom: + 0 exactly one candidate is valid + 1 both candidates are valid + 2 neither candidate is valid (includes coplanar/parallel case) """ px, py, pz, ngx, ngy, ngz = _accux_gca_scalar( w00, w01, w02, w10, w11, w12, v00, v01, v02, v10, v11, v12 ) pos_fin = ( - int(math.isfinite(px)) * int(math.isfinite(py)) * int(math.isfinite(pz)) + int(math.isfinite(px)) + * int(math.isfinite(py)) + * int(math.isfinite(pz)) ) neg_fin = ( - int(math.isfinite(ngx)) * int(math.isfinite(ngy)) * int(math.isfinite(ngz)) + int(math.isfinite(ngx)) + * int(math.isfinite(ngy)) + * int(math.isfinite(ngz)) ) pos_on_a = pos_fin * _on_minor_arc_xyz(px, py, pz, w00, w01, w02, w10, w11, w12) pos_on_b = pos_fin * _on_minor_arc_xyz(px, py, pz, v00, v01, v02, v10, v11, v12) From 309d8051d14bc6d1a7bf425e414de12aed60c6ab Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 19 Aug 2026 22:03:06 +0000 Subject: [PATCH 4/6] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- uxarray/grid/intersections.py | 14 +++----------- 1 file changed, 3 insertions(+), 11 deletions(-) diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py index 1a8818bd4..3f4201f80 100644 --- a/uxarray/grid/intersections.py +++ b/uxarray/grid/intersections.py @@ -314,9 +314,7 @@ def _accux_gca(w0, w1, v0, v1): @njit(cache=True, inline="always", error_model="numpy") -def _accux_gca_scalar( - w00, w01, w02, w10, w11, w12, v00, v01, v02, v10, v11, v12 -): +def _accux_gca_scalar(w00, w01, w02, w10, w11, w12, v00, v01, v02, v10, v11, v12): """Scalar-argument form of :func:`_accux_gca`: returns the six components. Compute the candidate intersection points of two great-circle arcs. @@ -442,15 +440,9 @@ def _try_gca_gca_intersection_scalar( w00, w01, w02, w10, w11, w12, v00, v01, v02, v10, v11, v12 ) - pos_fin = ( - int(math.isfinite(px)) - * int(math.isfinite(py)) - * int(math.isfinite(pz)) - ) + pos_fin = int(math.isfinite(px)) * int(math.isfinite(py)) * int(math.isfinite(pz)) neg_fin = ( - int(math.isfinite(ngx)) - * int(math.isfinite(ngy)) - * int(math.isfinite(ngz)) + int(math.isfinite(ngx)) * int(math.isfinite(ngy)) * int(math.isfinite(ngz)) ) pos_on_a = pos_fin * _on_minor_arc_xyz(px, py, pz, w00, w01, w02, w10, w11, w12) pos_on_b = pos_fin * _on_minor_arc_xyz(px, py, pz, v00, v01, v02, v10, v11, v12) From 97c044be362c2efb1ca988bdaa36188e265f9556 Mon Sep 17 00:00:00 2001 From: cmdupuis3 Date: Mon, 24 Aug 2026 16:57:27 -0500 Subject: [PATCH 5/6] scalarized gca intersections deslop --- benchmarks/geometry_kernels.py | 45 ++----- uxarray/grid/intersections.py | 210 +++++++++++++++------------------ 2 files changed, 108 insertions(+), 147 deletions(-) diff --git a/benchmarks/geometry_kernels.py b/benchmarks/geometry_kernels.py index 5750bed46..9744704fa 100644 --- a/benchmarks/geometry_kernels.py +++ b/benchmarks/geometry_kernels.py @@ -14,7 +14,6 @@ """ import numpy as np -from asv_runner.benchmarks.mark import skip_benchmark_if def _unit(v): @@ -143,54 +142,32 @@ def time_on_minor_arc(self): self.on_minor_arc(_V0, _W0, _W1) -# ``gca_gca_intersection`` (Layer 3) calls scalar-argument L1/L2 kernels -# (``_accux_gca_scalar``, ``_try_gca_gca_intersection_scalar``) rather than the -# array-returning ``_accux_gca``/``_try_gca_gca_intersection`` that predate them -# -- the array forms each allocated a ``(3,)`` (or ``(2,3)``) result per call, -# which no longer happens on the hot path. Benchmarking the array forms would -# time functions the dispatcher does not call, so they are skipped -- not -# measured -- on commits before the scalar kernels existed. -try: - from uxarray.grid.intersections import ( # noqa: F401 - _accux_gca_scalar, - _try_gca_gca_intersection_scalar, - ) - - _HAS_SCALAR_GCA_KERNELS = True -except ImportError: - _HAS_SCALAR_GCA_KERNELS = False - - class GCAGCAIntersection: """Benchmark all three layers of the GCA-GCA intersection stack.""" def setup(self): - from uxarray.grid.intersections import gca_gca_intersection + from uxarray.grid.intersections import ( + _accux_gca_scalar, + _try_gca_gca_intersection_scalar, + gca_gca_intersection + ) self.gca_gca_intersection = gca_gca_intersection self.gca_a = np.stack([_W0, _W1]) self.gca_b = np.stack([_V0, _V1]) gca_gca_intersection(self.gca_a, self.gca_b) - if _HAS_SCALAR_GCA_KERNELS: - from uxarray.grid.intersections import ( - _accux_gca_scalar, - _try_gca_gca_intersection_scalar, - ) - - self._accux_gca_scalar = _accux_gca_scalar - self._try_gca_gca_intersection_scalar = _try_gca_gca_intersection_scalar - _accux_gca_scalar(*_W0, *_W1, *_V0, *_V1) - _try_gca_gca_intersection_scalar(*_W0, *_W1, *_V0, *_V1) + self._accux_gca_scalar = _accux_gca_scalar + self._try_gca_gca_intersection_scalar = _try_gca_gca_intersection_scalar + _accux_gca_scalar(*_W0, *_W1, *_V0, *_V1) + _try_gca_gca_intersection_scalar(*_W0, *_W1, *_V0, *_V1) - @skip_benchmark_if(not _HAS_SCALAR_GCA_KERNELS) def time_accux_gca_kernel(self): - """Layer 1: pure numerical kernel (scalar form; allocation-free).""" + """Layer 1: pure numerical kernel.""" self._accux_gca_scalar(*_W0, *_W1, *_V0, *_V1) - @skip_benchmark_if(not _HAS_SCALAR_GCA_KERNELS) def time_try_gca_gca_intersection(self): - """Layer 2: batch/status layer (scalar form; allocation-free).""" + """Layer 2: batch/status layer.""" self._try_gca_gca_intersection_scalar(*_W0, *_W1, *_V0, *_V1) def time_gca_gca_intersection(self): diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py index 3f4201f80..e7debdc6e 100644 --- a/uxarray/grid/intersections.py +++ b/uxarray/grid/intersections.py @@ -251,68 +251,6 @@ def faces_within_lat_bounds(lats, face_bounds_lat): return _flatnonzero(within_bounds) -@njit(cache=True, inline="always", error_model="numpy") -def _accux_gca(w0, w1, v0, v1): - """Compute the candidate intersection points of two great-circle arcs. - - Pure numerical kernel (mirrors AccuSphGeom ``accux_gca``). - - Computes the two antipodal candidate intersection points of the great-circle - arcs w0-w1 and v0-v1. No branching, no validity filtering. - - Parameters - ---------- - w0, w1 : np.ndarray, shape (3,) - Cartesian endpoints of the first arc. - v0, v1 : np.ndarray, shape (3,) - Cartesian endpoints of the second arc. - - Returns - ------- - pos, neg : np.ndarray, shape (3,) - Two antipodal candidate unit vectors. - """ - n1x_hi, n1y_hi, n1z_hi, n1x_lo, n1y_lo, n1z_lo = accucross( - w0[0], w0[1], w0[2], w1[0], w1[1], w1[2] - ) - n2x_hi, n2y_hi, n2z_hi, n2x_lo, n2y_lo, n2z_lo = accucross( - v0[0], v0[1], v0[2], v1[0], v1[1], v1[2] - ) - vx_hi, vy_hi, vz_hi, vx_lo, vy_lo, vz_lo = accucross_pair( - n1x_hi, - n1y_hi, - n1z_hi, - n1x_lo, - n1y_lo, - n1z_lo, - n2x_hi, - n2y_hi, - n2z_hi, - n2x_lo, - n2y_lo, - n2z_lo, - ) - vx = vx_hi + vx_lo - vy = vy_hi + vy_lo - vz = vz_hi + vz_lo - # Compensated norm: sum_of_squares_c over the (hi, lo) vector, then acc_sqrt_re - # folding the low part into the root, matching AccuSphGeom accux_gca. n = root.hi. - sum_hi, sum_lo = _sum_of_squares_c((vx_hi, vy_hi, vz_hi), (vx_lo, vy_lo, vz_lo)) - vn, _ = acc_sqrt_re(sum_hi, sum_lo) - # vn==0 (coplanar arcs) yields inf via IEEE division under error_model="numpy", - # so pos/neg become non-finite and the status layer masks them out. Branch-free. - inv = 1.0 / vn - pos = np.empty(3) - pos[0] = vx * inv - pos[1] = vy * inv - pos[2] = vz * inv - neg = np.empty(3) - neg[0] = -pos[0] - neg[1] = -pos[1] - neg[2] = -pos[2] - return pos, neg - - @njit(cache=True, inline="always", error_model="numpy") def _accux_gca_scalar(w00, w01, w02, w10, w11, w12, v00, v01, v02, v10, v11, v12): """Scalar-argument form of :func:`_accux_gca`: returns the six components. @@ -359,6 +297,8 @@ def _accux_gca_scalar(w00, w01, w02, w10, w11, w12, v00, v01, v02, v10, v11, v12 vx = vx_hi + vx_lo vy = vy_hi + vy_lo vz = vz_hi + vz_lo + # Compensated norm: sum_of_squares_c over the (hi, lo) vector, then acc_sqrt_re + # folding the low part into the root, matching AccuSphGeom accux_gca. n = root.hi. sum_hi, sum_lo = _sum_of_squares_c((vx_hi, vy_hi, vz_hi), (vx_lo, vy_lo, vz_lo)) vn, _ = acc_sqrt_re(sum_hi, sum_lo) # vn==0 (coplanar arcs) yields inf via IEEE division under error_model="numpy", @@ -370,52 +310,50 @@ def _accux_gca_scalar(w00, w01, w02, w10, w11, w12, v00, v01, v02, v10, v11, v12 return pos_x, pos_y, pos_z, -pos_x, -pos_y, -pos_z -@njit(cache=True, error_model="numpy") -def _try_gca_gca_intersection(w0, w1, v0, v1): - """Select the valid great-circle intersection and report a status code. +@njit(cache=True, inline="always", error_model="numpy") +def _accux_gca(w0, w1, v0, v1): + """Compute the candidate intersection points of two great-circle arcs. - Batch/status layer (mirrors AccuSphGeom ``try_gca_gca_intersection``). + Pure numerical kernel (mirrors AccuSphGeom ``accux_gca``). - Calls the pure numerical kernel, applies integer mask arithmetic to determine - validity, selects the output point without if/else branching in the hot path. + Computes the two antipodal candidate intersection points of the great-circle + arcs w0-w1 and v0-v1. No branching, no validity filtering. - Status codes mirror AccuSphGeom: - 0 exactly one candidate is valid - 1 both candidates are valid - 2 neither candidate is valid (includes coplanar/parallel case) - """ - pos, neg = _accux_gca(w0, w1, v0, v1) + Parameters + ---------- + w0, w1 : np.ndarray, shape (3,) + Cartesian endpoints of the first arc. + v0, v1 : np.ndarray, shape (3,) + Cartesian endpoints of the second arc. - pos_fin = ( - int(math.isfinite(pos[0])) - * int(math.isfinite(pos[1])) - * int(math.isfinite(pos[2])) - ) - neg_fin = ( - int(math.isfinite(neg[0])) - * int(math.isfinite(neg[1])) - * int(math.isfinite(neg[2])) + Returns + ------- + pos, neg : np.ndarray, shape (3,) + Two antipodal candidate unit vectors. + """ + pos_x, pos_y, pos_z, neg_x, neg_y, neg_z = _accux_gca_scalar( + w0[0], + w0[1], + w0[2], + w1[0], + w1[1], + w1[2], + v0[0], + v0[1], + v0[2], + v1[0], + v1[1], + v1[2], ) - pos_on_a = pos_fin * on_minor_arc(pos, w0, w1) - pos_on_b = pos_fin * on_minor_arc(pos, v0, v1) - neg_on_a = neg_fin * on_minor_arc(neg, w0, w1) - neg_on_b = neg_fin * on_minor_arc(neg, v0, v1) - - pos_valid = pos_fin * pos_on_a * pos_on_b - neg_valid = neg_fin * neg_on_a * neg_on_b - - pos_mask = pos_valid * (1 - neg_valid) - neg_mask = neg_valid * (1 - pos_valid) - - point = np.empty(3) - point[0] = pos_mask * pos[0] + neg_mask * neg[0] - point[1] = pos_mask * pos[1] + neg_mask * neg[1] - point[2] = pos_mask * pos[2] + neg_mask * neg[2] - - both = pos_valid * neg_valid - none = (1 - pos_valid) * (1 - neg_valid) - status = both + none * 2 - return point, status, pos, neg + pos = np.empty(3) + pos[0] = pos_x + pos[1] = pos_y + pos[2] = pos_z + neg = np.empty(3) + neg[0] = neg_x + neg[1] = neg_y + neg[2] = neg_z + return pos, neg @njit(cache=True, inline="always", error_model="numpy") @@ -465,6 +403,61 @@ def _try_gca_gca_intersection_scalar( return point_x, point_y, point_z, status, px, py, pz, ngx, ngy, ngz +@njit(cache=True, error_model="numpy") +def _try_gca_gca_intersection(w0, w1, v0, v1): + """Select the valid great-circle intersection and report a status code. + + Batch/status layer (mirrors AccuSphGeom ``try_gca_gca_intersection``). + + Calls the pure numerical kernel, applies integer mask arithmetic to determine + validity, selects the output point without if/else branching in the hot path. + + Status codes mirror AccuSphGeom: + 0 exactly one candidate is valid + 1 both candidates are valid + 2 neither candidate is valid (includes coplanar/parallel case) + """ + ( + point_x, + point_y, + point_z, + status, + pos_x, + pos_y, + pos_z, + neg_x, + neg_y, + neg_z, + ) = _try_gca_gca_intersection_scalar( + w0[0], + w0[1], + w0[2], + w1[0], + w1[1], + w1[2], + v0[0], + v0[1], + v0[2], + v1[0], + v1[1], + v1[2], + ) + + point = np.empty(3) + point[0] = point_x + point[1] = point_y + point[2] = point_z + pos = np.empty(3) + pos[0] = pos_x + pos[1] = pos_y + pos[2] = pos_z + neg = np.empty(3) + neg[0] = neg_x + neg[1] = neg_y + neg[2] = neg_z + return point, status, pos, neg + + @njit(cache=True, error_model="numpy") def gca_gca_intersection(gca_a_xyz, gca_b_xyz): """Return the intersection points of two great-circle arcs. @@ -500,12 +493,7 @@ def gca_gca_intersection(gca_a_xyz, gca_b_xyz): if gca_a_xyz.shape[1] != 3 or gca_b_xyz.shape[1] != 3: raise DimensionError("The two GCAs must be in the cartesian [x, y, z] format") - # Unpack to scalars and run the allocation-free scalar chain. The array - # forms (_accux_gca, _try_gca_gca_intersection) allocated four (3,)/(2,3) - # arrays per call -- pos, neg, point, res -- which dominated this function - # in the face-bounds path, where it runs once per edge of every face. Only - # the (2, 3) result array remains, because the public return type is an - # array. The arithmetic is unchanged, so results are bit-for-bit identical. + # Unpack to scalars and run the allocation-free scalar chain. w00 = gca_a_xyz[0, 0] w01 = gca_a_xyz[0, 1] w02 = gca_a_xyz[0, 2] @@ -536,10 +524,6 @@ def gca_gca_intersection(gca_a_xyz, gca_b_xyz): res = np.empty((2, 3)) count = 0 - # The branch on ``status`` is load-bearing, not stylistic: when no candidate - # is valid both are non-finite, and the masked ``point_*`` above is nan. A - # mask-arithmetic rewrite of this dispatch would propagate that nan into the - # coplanar result (0.0 * nan == nan) instead of discarding it. if status == 0: res[0, 0] = point_x res[0, 1] = point_y From 0b83f9a92fefa2b634684d988efdfa583c6d8118 Mon Sep 17 00:00:00 2001 From: cmdupuis3 Date: Mon, 24 Aug 2026 19:06:36 -0500 Subject: [PATCH 6/6] Remove vector-valued routines --- benchmarks/geometry_kernels.py | 16 ++--- uxarray/grid/intersections.py | 117 ++------------------------------- 2 files changed, 14 insertions(+), 119 deletions(-) diff --git a/benchmarks/geometry_kernels.py b/benchmarks/geometry_kernels.py index 9744704fa..3892c0911 100644 --- a/benchmarks/geometry_kernels.py +++ b/benchmarks/geometry_kernels.py @@ -147,8 +147,8 @@ class GCAGCAIntersection: def setup(self): from uxarray.grid.intersections import ( - _accux_gca_scalar, - _try_gca_gca_intersection_scalar, + _accux_gca, + _try_gca_gca_intersection, gca_gca_intersection ) @@ -157,18 +157,18 @@ def setup(self): self.gca_b = np.stack([_V0, _V1]) gca_gca_intersection(self.gca_a, self.gca_b) - self._accux_gca_scalar = _accux_gca_scalar - self._try_gca_gca_intersection_scalar = _try_gca_gca_intersection_scalar - _accux_gca_scalar(*_W0, *_W1, *_V0, *_V1) - _try_gca_gca_intersection_scalar(*_W0, *_W1, *_V0, *_V1) + self._accux_gca = _accux_gca + self._try_gca_gca_intersection = _try_gca_gca_intersection + _accux_gca(*_W0, *_W1, *_V0, *_V1) + _try_gca_gca_intersection(*_W0, *_W1, *_V0, *_V1) def time_accux_gca_kernel(self): """Layer 1: pure numerical kernel.""" - self._accux_gca_scalar(*_W0, *_W1, *_V0, *_V1) + self._accux_gca(*_W0, *_W1, *_V0, *_V1) def time_try_gca_gca_intersection(self): """Layer 2: batch/status layer.""" - self._try_gca_gca_intersection_scalar(*_W0, *_W1, *_V0, *_V1) + self._try_gca_gca_intersection(*_W0, *_W1, *_V0, *_V1) def time_gca_gca_intersection(self): """Layer 3: dispatcher (full public API).""" diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py index e7debdc6e..d858eca06 100644 --- a/uxarray/grid/intersections.py +++ b/uxarray/grid/intersections.py @@ -252,10 +252,8 @@ def faces_within_lat_bounds(lats, face_bounds_lat): @njit(cache=True, inline="always", error_model="numpy") -def _accux_gca_scalar(w00, w01, w02, w10, w11, w12, v00, v01, v02, v10, v11, v12): - """Scalar-argument form of :func:`_accux_gca`: returns the six components. - - Compute the candidate intersection points of two great-circle arcs. +def _accux_gca(w00, w01, w02, w10, w11, w12, v00, v01, v02, v10, v11, v12): + """Compute the candidate intersection points of two great-circle arcs. Pure numerical kernel (mirrors AccuSphGeom ``accux_gca``). @@ -311,58 +309,10 @@ def _accux_gca_scalar(w00, w01, w02, w10, w11, w12, v00, v01, v02, v10, v11, v12 @njit(cache=True, inline="always", error_model="numpy") -def _accux_gca(w0, w1, v0, v1): - """Compute the candidate intersection points of two great-circle arcs. - - Pure numerical kernel (mirrors AccuSphGeom ``accux_gca``). - - Computes the two antipodal candidate intersection points of the great-circle - arcs w0-w1 and v0-v1. No branching, no validity filtering. - - Parameters - ---------- - w0, w1 : np.ndarray, shape (3,) - Cartesian endpoints of the first arc. - v0, v1 : np.ndarray, shape (3,) - Cartesian endpoints of the second arc. - - Returns - ------- - pos, neg : np.ndarray, shape (3,) - Two antipodal candidate unit vectors. - """ - pos_x, pos_y, pos_z, neg_x, neg_y, neg_z = _accux_gca_scalar( - w0[0], - w0[1], - w0[2], - w1[0], - w1[1], - w1[2], - v0[0], - v0[1], - v0[2], - v1[0], - v1[1], - v1[2], - ) - pos = np.empty(3) - pos[0] = pos_x - pos[1] = pos_y - pos[2] = pos_z - neg = np.empty(3) - neg[0] = neg_x - neg[1] = neg_y - neg[2] = neg_z - return pos, neg - - -@njit(cache=True, inline="always", error_model="numpy") -def _try_gca_gca_intersection_scalar( +def _try_gca_gca_intersection( w00, w01, w02, w10, w11, w12, v00, v01, v02, v10, v11, v12 ): - """Scalar-argument form of :func:`_try_gca_gca_intersection`. - - Select the valid great-circle intersection and report a status code. + """Select the valid great-circle intersection and report a status code. Batch/status layer (mirrors AccuSphGeom ``try_gca_gca_intersection``). @@ -374,7 +324,7 @@ def _try_gca_gca_intersection_scalar( 1 both candidates are valid 2 neither candidate is valid (includes coplanar/parallel case) """ - px, py, pz, ngx, ngy, ngz = _accux_gca_scalar( + px, py, pz, ngx, ngy, ngz = _accux_gca( w00, w01, w02, w10, w11, w12, v00, v01, v02, v10, v11, v12 ) @@ -403,61 +353,6 @@ def _try_gca_gca_intersection_scalar( return point_x, point_y, point_z, status, px, py, pz, ngx, ngy, ngz -@njit(cache=True, error_model="numpy") -def _try_gca_gca_intersection(w0, w1, v0, v1): - """Select the valid great-circle intersection and report a status code. - - Batch/status layer (mirrors AccuSphGeom ``try_gca_gca_intersection``). - - Calls the pure numerical kernel, applies integer mask arithmetic to determine - validity, selects the output point without if/else branching in the hot path. - - Status codes mirror AccuSphGeom: - 0 exactly one candidate is valid - 1 both candidates are valid - 2 neither candidate is valid (includes coplanar/parallel case) - """ - ( - point_x, - point_y, - point_z, - status, - pos_x, - pos_y, - pos_z, - neg_x, - neg_y, - neg_z, - ) = _try_gca_gca_intersection_scalar( - w0[0], - w0[1], - w0[2], - w1[0], - w1[1], - w1[2], - v0[0], - v0[1], - v0[2], - v1[0], - v1[1], - v1[2], - ) - - point = np.empty(3) - point[0] = point_x - point[1] = point_y - point[2] = point_z - pos = np.empty(3) - pos[0] = pos_x - pos[1] = pos_y - pos[2] = pos_z - neg = np.empty(3) - neg[0] = neg_x - neg[1] = neg_y - neg[2] = neg_z - return point, status, pos, neg - - @njit(cache=True, error_model="numpy") def gca_gca_intersection(gca_a_xyz, gca_b_xyz): """Return the intersection points of two great-circle arcs. @@ -518,7 +413,7 @@ def gca_gca_intersection(gca_a_xyz, gca_b_xyz): neg_x, neg_y, neg_z, - ) = _try_gca_gca_intersection_scalar( + ) = _try_gca_gca_intersection( w00, w01, w02, w10, w11, w12, v00, v01, v02, v10, v11, v12 )