diff --git a/.github/workflows/asv-benchmarking-comment.yml b/.github/workflows/asv-benchmarking-comment.yml index 0b6befbd4..a6852600f 100644 --- a/.github/workflows/asv-benchmarking-comment.yml +++ b/.github/workflows/asv-benchmarking-comment.yml @@ -24,7 +24,7 @@ jobs: if: ${{ github.event.workflow_run.conclusion == 'success' }} steps: - name: Download benchmark results artifact - uses: actions/download-artifact@v7 + uses: actions/download-artifact@v8 with: name: asv-benchmark-results-Linux path: asv-results diff --git a/.github/workflows/yac-optional.yml b/.github/workflows/yac-optional.yml index 68b166bbb..ccfa606db 100644 --- a/.github/workflows/yac-optional.yml +++ b/.github/workflows/yac-optional.yml @@ -10,13 +10,13 @@ on: jobs: yac-optional: - name: YAC core v3.14.0_p1 (Ubuntu) + name: YAC core v3.18.0 (Ubuntu) runs-on: ubuntu-latest defaults: run: shell: bash -l {0} env: - YAC_VERSION: v3.14.0_p1 + YAC_VERSION: v3.18.0 YAXT_VERSION: v0.11.5.1 MPIEXEC: /usr/bin/mpirun MPIRUN: /usr/bin/mpirun @@ -36,7 +36,7 @@ jobs: with: activate-environment: uxarray_build channel-priority: strict - python-version: "3.11" + python-version: "3.14" channels: conda-forge environment-file: ci/environment.yml miniforge-variant: Miniforge3 @@ -65,8 +65,9 @@ jobs: mpif90 --version - name: Install Python build dependencies run: | - python -m pip install --upgrade pip - python -m pip install cython wheel + # Cython and wheel are needed to build YAC's Python bindings from source. + conda install --yes -c conda-forge "cython>=3.1" wheel + python -m cython --version - name: Build and install YAXT run: | set -euxo pipefail diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index edcf0ed53..838b780f7 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -14,7 +14,7 @@ repos: - repo: https://github.com/astral-sh/ruff-pre-commit # Ruff version. - rev: v0.15.20 + rev: v0.15.21 hooks: - id: ruff name: lint with ruff diff --git a/benchmarks/geometry_samebody.py b/benchmarks/geometry_samebody.py index 5f9e96d46..7712d1c7a 100644 --- a/benchmarks/geometry_samebody.py +++ b/benchmarks/geometry_samebody.py @@ -201,21 +201,28 @@ def _unit(v): return v / np.linalg.norm(v) -def _make_cases(n, seed): - """Random great-circle arcs paired with a latitude their arc actually crosses.""" +def _make_cases(n, seed, frac_cross=0.7): + rng = np.random.default_rng(seed) cases = [] while len(cases) < n: - a = _unit(rng.standard_normal(3)) - b = _unit(rng.standard_normal(3)) - if abs(np.dot(a, b)) > 0.999: # near-degenerate arc, skip - continue - # pick a latitude strictly between the two endpoints' z so an - # intersection is likely to exist - zlo, zhi = sorted((a[2], b[2])) - if zhi - zlo < 1e-6: - continue - const_z = zlo + (zhi - zlo) * rng.uniform(0.2, 0.8) + mid = _unit(rng.standard_normal(3)) + tan = _unit(np.cross(mid, _unit(rng.standard_normal(3)))) + half = 0.5 * math.radians(10.0 ** rng.uniform(math.log10(0.2), math.log10(8.0))) + ca, sa = math.cos(half), math.sin(half) + a = _unit(mid * ca - tan * sa) + b = _unit(mid * ca + tan * sa) + zlo, zhi = (a[2], b[2]) if a[2] <= b[2] else (b[2], a[2]) + if rng.random() < frac_cross and zhi - zlo > 1e-9: + const_z = zlo + (zhi - zlo) * rng.uniform(0.15, 0.85) + else: + # a latitude the arc does not span -> empty-return path + room_above = 0.999 - zhi + room_below = zlo + 0.999 + if room_above >= room_below: + const_z = zhi + min(room_above, rng.uniform(0.02, 0.3)) + else: + const_z = zlo - min(room_below, rng.uniform(0.02, 0.3)) cases.append((np.stack([a, b]), float(const_z))) return cases @@ -324,10 +331,8 @@ def main(): if np.isfinite(d): max_out_diff = max(max_out_diff, d) - # ---- timing: replicate the 200 cases into a large batch (~200k points) ---- - reps = 1000 - big = base_cases * reps - A, B, Z, gcas = _pack(big) + # ---- timing: a large batch of DISTINCT cases (no replication) ---- + A, B, Z, gcas = _pack(_make_cases(100_000, seed=20251105)) n = A.shape[0] t_direct = _time_batch(_batch_fp64_kernel, (A, B, Z)) @@ -341,7 +346,7 @@ def main(): print("Same-body FP64-vs-AccuX diagnostic — GCA/ConstLat (PR #1513)") print("=" * 70) print(f"baseline cases: {len(base_cases)} with-result: {n_with_result}") - print(f"timing batch : {n} points ({reps}x replication), best of 7") + print(f"timing batch : {n} distinct points, best of 7") print() print("CORRECTNESS (same-body FP64 dispatcher vs real AccuX dispatcher)") print(f" status mismatches : {status_mismatch}") @@ -391,7 +396,7 @@ class SameBodyConstLat: """asv: same-body FP64 vs real AccuX at kernel (L1) and dispatch (L3) levels.""" def setup(self): - cases = _make_cases(200, seed=20251104) * 100 + cases = _make_cases(20_000, seed=20251104) self.A, self.B, self.Z, self.gcas = _pack(cases) _batch_fp64_kernel(self.A, self.B, self.Z) _batch_accux_kernel(self.A, self.B, self.Z) diff --git a/benchmarks/geometry_samebody_gcagca.py b/benchmarks/geometry_samebody_gcagca.py new file mode 100644 index 000000000..348f55859 --- /dev/null +++ b/benchmarks/geometry_samebody_gcagca.py @@ -0,0 +1,338 @@ +"""Same-body FP64-vs-AccuX diagnostic for the GCA x GCA path. + +Companion to ``geometry_samebody.py``, which only covers the GCA/constant-latitude +stack. T + +gca_gca is the kernel invoked per edge from the njit point-in-polygon crossing +loop (``uxarray/grid/geometry.py`` ``_check_intersection``), so its allocation +profile is on the hot path of ``get_faces_containing_point`` and pole detection. + +Factoring (identical to ``geometry_samebody.py``):: + + real AccuX time - same-body FP64 time == EFT math only + same-body FP64 - direct FP64 kernel == L2/L3 plumbing only + +""" + +import math +import time + +import numpy as np +from numba import njit + +from uxarray.grid.arcs import on_minor_arc +from uxarray.grid.intersections import _accux_gca, gca_gca_intersection + + +@njit(cache=True, inline="always") +def _fp64_gca(w0, w1, v0, v1): + """ + L1 (FP64 body) -- plain double-precision cross-product triple, the direct + analogue of _accux_gca (intersections.py) with accucross/accucross_pair + replaced by naive FP64 cross products. Same allocation shape (two np.empty(3)) + so the twin's allocation profile matches the real kernel exactly. + """ + + n1x = w0[1] * w1[2] - w0[2] * w1[1] + n1y = w0[2] * w1[0] - w0[0] * w1[2] + n1z = w0[0] * w1[1] - w0[1] * w1[0] + n2x = v0[1] * v1[2] - v0[2] * v1[1] + n2y = v0[2] * v1[0] - v0[0] * v1[2] + n2z = v0[0] * v1[1] - v0[1] * v1[0] + + vx = n1y * n2z - n1z * n2y + vy = n1z * n2x - n1x * n2z + vz = n1x * n2y - n1y * n2x + + vn = math.sqrt(vx * vx + vy * vy + vz * vz) + inv = 1.0 / vn if vn != 0.0 else np.inf + 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) +def _fp64_try_gca_gca_intersection(w0, w1, v0, v1): + """ + L2 (FP64 body) -- byte-for-byte identical logic to _try_gca_gca_intersection + (intersections.py), only the L1 call differs. + """ + pos, neg = _fp64_gca(w0, w1, v0, v1) + + pos_fin = ( + 1 + if math.isfinite(pos[0]) and math.isfinite(pos[1]) and math.isfinite(pos[2]) + else 0 + ) + neg_fin = ( + 1 + if math.isfinite(neg[0]) and math.isfinite(neg[1]) and math.isfinite(neg[2]) + else 0 + ) + pos_on_a = 1 if (pos_fin and on_minor_arc(pos, w0, w1)) else 0 + pos_on_b = 1 if (pos_fin and on_minor_arc(pos, v0, v1)) else 0 + neg_on_a = 1 if (neg_fin and on_minor_arc(neg, w0, w1)) else 0 + neg_on_b = 1 if (neg_fin and on_minor_arc(neg, v0, v1)) else 0 + + 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 + + +@njit(cache=True) +def _fp64_gca_gca_intersection(gca_a_xyz, gca_b_xyz): + """ + L3 (FP64 body) -- identical dispatcher to gca_gca_intersection + (intersections.py), same np.empty((2, 3)) + res[:count] slice profile. + """ + if gca_a_xyz.shape[1] != 3 or gca_b_xyz.shape[1] != 3: + raise ValueError("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 = _fp64_try_gca_gca_intersection(w0, w1, v0, v1) + + res = np.empty((2, 3)) + count = 0 + if status == 0: + res[0, 0] = point[0] + res[0, 1] = point[1] + res[0, 2] = point[2] + 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] + count = 2 + else: + if on_minor_arc(v0, w0, w1): + res[count, 0] = v0[0] + res[count, 1] = v0[1] + res[count, 2] = v0[2] + count += 1 + if on_minor_arc(v1, w0, w1): + res[count, 0] = v1[0] + res[count, 1] = v1[1] + res[count, 2] = v1[2] + count += 1 + return res[:count] + + +# --------------------------------------------------------------------------- +# Case generation -- a LARGE set of DISTINCT, SHORT arcs (matching real grid +# edges, ~0.2-8 deg / median ~1.5 deg) with a controlled intersect / no-intersect +# mix. Short arcs put a x b in the near-parallel cancellation regime the EFT +# targets; the mix exercises both the status 0/1 (allocating) and status 2 +# (empty-return) dispatcher branches. +# --------------------------------------------------------------------------- + + +def _unit(v): + return v / np.linalg.norm(v) + + +def _short_arc(rng, mid=None): + """A short great-circle arc (~0.2-8 deg, median ~1.5 deg) around ``mid`` (a + random point if not given), matching real unstructured-grid edge lengths.""" + if mid is None: + mid = _unit(rng.standard_normal(3)) + tan = _unit(np.cross(mid, _unit(rng.standard_normal(3)))) + half = 0.5 * math.radians(10.0 ** rng.uniform(math.log10(0.2), math.log10(8.0))) + ca, sa = math.cos(half), math.sin(half) + return _unit(mid * ca - tan * sa), _unit(mid * ca + tan * sa) + + +def _make_gca_cases(n, seed, frac_intersect=0.6): + rng = np.random.default_rng(seed) + cases = [] + while len(cases) < n: + if rng.random() < frac_intersect: + # two short arcs sharing a midpoint p -> p is the midpoint of both + # minor arcs, so they intersect there (status 0/1). + p = _unit(rng.standard_normal(3)) + a0, a1 = _short_arc(rng, mid=p) + b0, b1 = _short_arc(rng, mid=p) + else: + # two independent short arcs -> their great circles cross off at + # least one minor arc, exercising the status 2 / empty-return branch. + a0, a1 = _short_arc(rng) + b0, b1 = _short_arc(rng) + cases.append((np.stack([a0, a1]), np.stack([b0, b1]))) + return cases + + +def _pack_gca(cases): + wa = np.ascontiguousarray([c[0][0] for c in cases]) + wb = np.ascontiguousarray([c[0][1] for c in cases]) + va = np.ascontiguousarray([c[1][0] for c in cases]) + vb = np.ascontiguousarray([c[1][1] for c in cases]) + ga = np.ascontiguousarray([c[0] for c in cases]) # (n, 2, 3) + gb = np.ascontiguousarray([c[1] for c in cases]) + return wa, wb, va, vb, ga, gb + + +# --------------------------------------------------------------------------- +# Batched in-kernel drivers -- the timing loop lives inside njit (mirrors +# geometry_samebody.py). Accumulate a scalar to defeat dead-code elimination. +# --------------------------------------------------------------------------- + + +@njit(cache=True) +def _batch_accux_gca_kernel(wa, wb, va, vb): + acc = 0.0 + for i in range(wa.shape[0]): + pos, neg = _accux_gca(wa[i], wb[i], va[i], vb[i]) + acc += pos[0] + pos[1] + pos[2] + return acc + + +@njit(cache=True) +def _batch_fp64_gca_kernel(wa, wb, va, vb): + acc = 0.0 + for i in range(wa.shape[0]): + pos, neg = _fp64_gca(wa[i], wb[i], va[i], vb[i]) + acc += pos[0] + pos[1] + pos[2] + return acc + + +@njit(cache=True) +def _batch_accux_gca_dispatch(ga, gb): + acc = 0.0 + for i in range(ga.shape[0]): + res = gca_gca_intersection(ga[i], gb[i]) + if res.shape[0] > 0: + acc += res[0, 0] + return acc + + +@njit(cache=True) +def _batch_fp64_gca_dispatch(ga, gb): + acc = 0.0 + for i in range(ga.shape[0]): + res = _fp64_gca_gca_intersection(ga[i], gb[i]) + if res.shape[0] > 0: + acc += res[0, 0] + return acc + + +def _time_batch(fn, args, repeat=7): + """Best-of-`repeat` wall-time for one batched call (compile excluded).""" + fn(*args) # warm / compile + best = math.inf + for _ in range(repeat): + t0 = time.perf_counter() + fn(*args) + best = min(best, time.perf_counter() - t0) + return best + + +def main(n_cases=100_000, seed=20251104): + """ + Standalone diagnostic driver (prints ns/edge-pair). Mirrors geometry_samebody.main. + """ + cases = _make_gca_cases(n_cases, seed=seed) + wa, wb, va, vb, ga, gb = _pack_gca(cases) + n = wa.shape[0] + + # correctness: FP64 twin vs real AccuX dispatcher (row count + max diff) + row_mismatch = 0 + max_out_diff = 0.0 + n_check = min(n, 5000) + for i in range(n_check): + r_fp = _fp64_gca_gca_intersection(ga[i], gb[i]) + r_ax = gca_gca_intersection(ga[i], gb[i]) + if r_fp.shape[0] != r_ax.shape[0]: + row_mismatch += 1 + elif r_fp.shape[0] > 0: + max_out_diff = max(max_out_diff, float(np.max(np.abs(r_fp - r_ax)))) + + t_fp64_k = _time_batch(_batch_fp64_gca_kernel, (wa, wb, va, vb)) + t_accux_k = _time_batch(_batch_accux_gca_kernel, (wa, wb, va, vb)) + t_fp64_d = _time_batch(_batch_fp64_gca_dispatch, (ga, gb)) + t_accux_d = _time_batch(_batch_accux_gca_dispatch, (ga, gb)) + + def ns(t): + return t / n * 1e9 + + print("=" * 70) + print("Same-body FP64-vs-AccuX diagnostic -- GCA x GCA") + print("=" * 70) + print(f"distinct cases : {n} (best of 7, in-kernel batch)") + print( + f"correctness : row mismatches {row_mismatch}/{n_check}" + f" max |diff| {max_out_diff:.3e}" + ) + print() + print("TIMING (ns per edge-pair)") + print(f" L1 FP64 kernel : {ns(t_fp64_k):8.2f} ns") + print(f" L1 AccuX kernel : {ns(t_accux_k):8.2f} ns") + print(f" L1+L2+L3 FP64 dispatch : {ns(t_fp64_d):8.2f} ns") + print(f" L1+L2+L3 AccuX dispatch : {ns(t_accux_d):8.2f} ns") + print() + print("DECOMPOSITION") + print(f" plumbing (FP64 body) : {ns(t_fp64_d - t_fp64_k):8.2f} ns/edge") + print(f" plumbing (AccuX body) : {ns(t_accux_d - t_accux_k):8.2f} ns/edge") + print( + f" EFT math (kernel) : {ns(t_accux_k - t_fp64_k):8.2f} ns/edge" + f" ({t_accux_k / t_fp64_k:.2f}x)" + ) + print( + f" EFT math (dispatch) : {ns(t_accux_d - t_fp64_d):8.2f} ns/edge" + f" ({t_accux_d / t_fp64_d:.2f}x)" + ) + print("=" * 70) + + +class SameBodyGcaGca: + """ + asv timing class (Numba warmed in setup, distinct cases) + same-body FP64 vs real AccuX gca_gca at kernel (L1) and dispatch (L3). + """ + + def setup(self): + cases = _make_gca_cases(100_000, seed=20251104) + self.wa, self.wb, self.va, self.vb, self.ga, self.gb = _pack_gca(cases) + _batch_fp64_gca_kernel(self.wa, self.wb, self.va, self.vb) + _batch_accux_gca_kernel(self.wa, self.wb, self.va, self.vb) + _batch_fp64_gca_dispatch(self.ga, self.gb) + _batch_accux_gca_dispatch(self.ga, self.gb) + + def time_fp64_kernel(self): + _batch_fp64_gca_kernel(self.wa, self.wb, self.va, self.vb) + + def time_accux_kernel(self): + _batch_accux_gca_kernel(self.wa, self.wb, self.va, self.vb) + + def time_fp64_dispatch(self): + _batch_fp64_gca_dispatch(self.ga, self.gb) + + def time_accux_dispatch(self): + _batch_accux_gca_dispatch(self.ga, self.gb) + + +if __name__ == "__main__": + main() diff --git a/ci/environment.yml b/ci/environment.yml index 2049e965e..44897b982 100644 --- a/ci/environment.yml +++ b/ci/environment.yml @@ -20,7 +20,6 @@ dependencies: - numba - numpy - pandas - - pathlib - pre_commit - polars - pyarrow diff --git a/pyproject.toml b/pyproject.toml index 3ee69940c..0596956c6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -9,6 +9,7 @@ classifiers=[ "Programming Language :: Python :: 3.11", "Programming Language :: Python :: 3.12", "Programming Language :: Python :: 3.13", + "Programming Language :: Python :: 3.14", ] dynamic = ["version"] @@ -30,7 +31,7 @@ dependencies = [ "matplotlib<3.11", # matplotlib 3.11.0 breaks some plots but <3.11 is fine; see issue #1542. "matplotlib-inline", "netcdf4", - "numba", + "numba>=0.63", # 0.63 is the first release supporting Python 3.14 (3.10 <= py < 3.15). See #1561. "numpy", "pandas", "pyarrow", @@ -50,7 +51,7 @@ dependencies = [ [project.optional-dependencies] complete = ["uxarray[dev]"] -dev = ['pathlib', 'pre_commit', 'pytest', 'pytest-cov', 'ruff', 'asv'] +dev = ['pre_commit', 'pytest', 'pytest-cov', 'ruff', 'asv'] [project.urls] Documentation = "https://uxarray.readthedocs.io/" diff --git a/test/grid/geometry/test_eft_accuracy.py b/test/grid/geometry/test_eft_accuracy.py new file mode 100644 index 000000000..58f831b4d --- /dev/null +++ b/test/grid/geometry/test_eft_accuracy.py @@ -0,0 +1,170 @@ +"""Accuracy regression tests for the compensated EFT primitives in +``uxarray.utils.computing`` — specifically the compensated sum-of-squares used by +the GCA / constant-latitude intersection kernel. + +An earlier revision computed the naive ``Σ h² + Σ l²`` instead of the +correct compensated squared norm ``Σ (h + l)² = Σ h² + 2·Σ h·l`` (matching +AccuSphGeom's ``numeric::sum_of_squares_c``. The correct compensated result +reaches ~1e-30 relative accuracy on well-conditioned inputs; the naive form +only reaches ~1e-15. +""" + +from fractions import Fraction + +import numpy as np +import pytest + +from uxarray.utils.computing import ( + _HAS_FMA, + _sum_of_squares_c, + _two_prod_veltkamp, + accucross, + two_prod, +) + +# Correct compensated result: ~1e-30 rel err. Naive Σh²+Σl²: ~1e-15. +_SUM_SQ_REL_TOL = 1e-24 + + +def _unit(v): + return v / np.linalg.norm(v) + + +def _make_normals(n, seed): + """Compensated cross products (hi, lo) of well-separated arc pairs. + + ``|dot(a, b)| < 0.999`` keeps ``|n|`` away from the noise floor so the + compensated squared norm is meaningful. (Extreme near-parallel arcs lose + accuracy for *any* algorithm and are out of scope for this primitive test.) + """ + rng = np.random.default_rng(seed) + hi = np.empty((n, 3)) + lo = np.empty((n, 3)) + i = 0 + while i < n: + a = _unit(rng.standard_normal(3)) + b = _unit(rng.standard_normal(3)) + if abs(np.dot(a, b)) > 0.999: + continue + xh, yh, zh, xl, yl, zl = accucross(a[0], a[1], a[2], b[0], b[1], b[2]) + hi[i] = (xh, yh, zh) + lo[i] = (xl, yl, zl) + i += 1 + return hi, lo + + +def _rel_err_vs_exact(res, hi_tup, lo_tup): + """|(res_hi + res_lo) − Σ (h + l)²| / Σ (h + l)², computed exactly.""" + true = sum((Fraction(h) + Fraction(l)) ** 2 for h, l in zip(hi_tup, lo_tup)) + if true == 0: + return 0.0 + return abs(float((Fraction(res[0]) + Fraction(res[1]) - true) / true)) + + +@pytest.fixture(scope="module") +def normals(): + return _make_normals(2000, seed=20260716) + + +@pytest.mark.parametrize("ncomp", [2, 3]) +def test_sum_of_squares_compensated_accuracy(normals, ncomp): + """The generic keeps the 2·Σh·l cross term for each tuple length used by the + const-lat kernel (N=2 -> denominator, N=3 -> |n|²). A naive Σh²+Σl² would + regress to ~1e-15 and fail this bound.""" + hi, lo = normals + worst = 0.0 + for i in range(hi.shape[0]): + h = tuple(float(x) for x in hi[i, :ncomp]) + lo_t = tuple(float(x) for x in lo[i, :ncomp]) + worst = max(worst, _rel_err_vs_exact(_sum_of_squares_c(h, lo_t), h, lo_t)) + assert worst < _SUM_SQ_REL_TOL, ( + f"N={ncomp}: _sum_of_squares_c max rel err {worst:.2e} ≥ " + f"{_SUM_SQ_REL_TOL:.0e} (naive Σh²+Σl² regressed the cross term?)" + ) + + +def test_sum_of_squares_keeps_cross_term(normals): + """Directly assert the result tracks the compensated ``Σh² + 2·Σh·l`` and + NOT the naive ``Σh² + Σl²``, on cases where the two formulas are exactly + distinguishable.""" + hi, lo = normals + checked = 0 + for i in range(hi.shape[0]): + h = (float(hi[i, 0]), float(hi[i, 1])) + lo_t = (float(lo[i, 0]), float(lo[i, 1])) + big_h = [Fraction(x) for x in h] + big_l = [Fraction(x) for x in lo_t] + correct = sum(big_h[k] * big_h[k] for k in range(2)) + 2 * sum( + big_h[k] * big_l[k] for k in range(2) + ) + naive = sum(big_h[k] * big_h[k] + big_l[k] * big_l[k] for k in range(2)) + # The two formulas differ by ~2·Σh·l ≈ 1e-16 relative — that IS the bug + # we test. Skip only cases where the gap falls below the compensated + # accuracy floor (so the two are genuinely indistinguishable there). + if correct == 0 or abs(float((naive - correct) / correct)) < 1e-20: + continue + res = _sum_of_squares_c(h, lo_t) + val = Fraction(res[0]) + Fraction(res[1]) + assert abs(val - correct) < abs(val - naive) + checked += 1 + assert checked > 100, "test did not exercise enough distinguishable cases" + + +def test_sum_of_squares_generic_any_n(): + """The single generic works for any tuple length — e.g. N=4, which no + hand-written ``_sum_sq_cN`` exists for — with the same compensated accuracy. + This is the point of one N-generic primitive instead of per-N helpers. + """ + rng = np.random.default_rng(7) + worst = 0.0 + for _ in range(2000): + h = tuple(float(x) for x in rng.standard_normal(4)) + # residual-scale low parts (~1 ulp of each high part) + lo_t = tuple(x * 1e-16 * float(rng.standard_normal()) for x in h) + worst = max(worst, _rel_err_vs_exact(_sum_of_squares_c(h, lo_t), h, lo_t)) + assert worst < _SUM_SQ_REL_TOL, ( + f"N=4: _sum_of_squares_c max rel err {worst:.2e} ≥ {_SUM_SQ_REL_TOL:.0e}" + ) + + +def _random_operands(rng, n): + for _ in range(n): + yield ( + float(rng.standard_normal() * rng.integers(1, 1 << 20)), + float(rng.standard_normal() * rng.integers(1, 1 << 20)), + ) + + +def test_two_prod_error_term_is_the_exact_residual(): + """``two_prod`` must return the EXACT rounding error of ``a * b``. + + This is the property ``computing._validate_fma`` exists to guarantee at + import time, asserted here directly against exact rational arithmetic so it + cannot silently lapse. A non-fused FMA lowering (or a non-compliant FMA) + yields ``e = 0.0``, which fails this outright. + + Note this is deliberately checked as ``p + e == a*b`` **exactly, over the + rationals** — NOT as the float expression ``p + e``, which rounds straight + back to ``p`` (since ``|e| <= ulp(p)/2``) and would make the assertion + vacuously true for any ``e`` whatsoever. + """ + rng = np.random.default_rng(20260716) + for a, b in _random_operands(rng, 20000): + p, e = two_prod(a, b) + assert Fraction(p) + Fraction(e) == Fraction(a) * Fraction(b), ( + f"two_prod({a!r}, {b!r}) = ({p!r}, {e!r}) is not an exact " + f"error-free transform" + ) + + +@pytest.mark.skipif(not _HAS_FMA, reason="no FMA path on this toolchain") +def test_two_prod_fma_matches_veltkamp_bit_for_bit(): + """The FMA and Veltkamp paths must agree bit-for-bit, so ``_HAS_FMA`` can + never change results. The exact residual is unique and representable, so + two correct implementations have no freedom to differ.""" + rng = np.random.default_rng(20260717) + for a, b in _random_operands(rng, 20000): + p, e = two_prod(a, b) + pv, ev = _two_prod_veltkamp(a, b) + assert np.float64(p).view(np.int64) == np.float64(pv).view(np.int64) + assert np.float64(e).view(np.int64) == np.float64(ev).view(np.int64) diff --git a/test/grid/grid/test_areas.py b/test/grid/grid/test_areas.py index 9f74d7d72..e3f1a54c6 100644 --- a/test/grid/grid/test_areas.py +++ b/test/grid/grid/test_areas.py @@ -30,6 +30,29 @@ def test_face_areas_calculate_total_face_area_triangle(mesh_constants): nt.assert_almost_equal(area_gaussian, mesh_constants['CORRECTED_TRI_AREA'], decimal=3) +def test_calculate_total_face_area_respects_quadrature_kwargs(mesh_constants): + """calculate_total_face_area must honor its quadrature_rule/order/ + latitude_adjusted_area kwargs instead of always returning the cached + default-parameter face areas.""" + verts = [ + [[0.02974582, -0.74469018, 0.66674712], + [0.1534193, -0.88744577, 0.43462917], + [0.18363692, -0.72230586, 0.66674712]] + ] + + grid_verts = ux.open_grid(verts, latlon=False) + + area_default = grid_verts.calculate_total_face_area( + quadrature_rule="triangular", order=4) + area_corrected = grid_verts.calculate_total_face_area( + quadrature_rule="gaussian", order=5, latitude_adjusted_area=True) + + # the corrected result must actually differ from the uncorrected one + assert not np.isclose(area_default, area_corrected) + nt.assert_almost_equal(area_default, mesh_constants['TRI_AREA'], decimal=6) + nt.assert_almost_equal(area_corrected, mesh_constants['CORRECTED_TRI_AREA'], decimal=6) + + def test_face_areas_compute_face_areas_geoflow_small(gridpath): """Checks if the GeoFlow Small can generate a face areas output.""" grid_geoflow = ux.open_grid(gridpath("ugrid", "geoflow-small", "grid.nc")) diff --git a/test/io/test_scrip.py b/test/io/test_scrip.py index 57b844802..4b30ddb0c 100644 --- a/test/io/test_scrip.py +++ b/test/io/test_scrip.py @@ -101,6 +101,25 @@ def test_open_multigrid_mask_active_value_default(gridpath): assert grids["atm"].n_face == expected_atm +def test_scrip_radians_units(gridpath): + """SCRIP files with coordinates in radians are converted to degrees on load.""" + # scrip_radians.nc has a 2-cell grid whose lat/lon are stored in radians. + # The expected degree values are: face_lon=[10, 20], face_lat=[30, 40]. + grid_file = gridpath("scrip", "scrip_radians", "scrip_radians_grid.nc") + grid = ux.open_grid(grid_file) + + expected_face_lon = np.array([10.0, 20.0]) + expected_face_lat = np.array([30.0, 40.0]) + # 7 unique nodes: the two cells share only one corner point (15, 35) + expected_node_lon = np.array([5., 5., 15., 15., 15., 25., 25.]) + expected_node_lat = np.array([25., 35., 25., 35., 45., 35., 45.]) + + nt.assert_allclose(np.sort(grid.face_lon.values), np.sort(expected_face_lon), atol=1e-10) + nt.assert_allclose(np.sort(grid.face_lat.values), np.sort(expected_face_lat), atol=1e-10) + nt.assert_allclose(np.sort(grid.node_lon.values), np.sort(expected_node_lon), atol=1e-10) + nt.assert_allclose(np.sort(grid.node_lat.values), np.sort(expected_node_lat), atol=1e-10) + + def test_open_multigrid_mask_active_value_per_grid_override(gridpath): """Per-grid override supports masks with different active values.""" grid_file = gridpath("scrip", "oasis", "grids.nc") diff --git a/test/meshfiles/scrip/scrip_radians/scrip_radians_grid.nc b/test/meshfiles/scrip/scrip_radians/scrip_radians_grid.nc new file mode 100644 index 000000000..437446ec7 Binary files /dev/null and b/test/meshfiles/scrip/scrip_radians/scrip_radians_grid.nc differ diff --git a/test/test_plot.py b/test/test_plot.py index 60c8a8d96..470cd63a7 100644 --- a/test/test_plot.py +++ b/test/test_plot.py @@ -249,6 +249,92 @@ def test_to_raster_pixel_ratio(gridpath, r1, r2): d = np.array(raster2.shape) - f * np.array(raster1.shape) assert (d >= 0).all() and (d <= f - 1).all() +def test_matplotlib_backend_restored_after_switch(monkeypatch): + """Regression test for #1537: switching HoloViews to matplotlib must restore + the Matplotlib backend captured beforehand. + + A bare pytest can't reproduce the Jupyter inline-hook flip, so we simulate + hv.extension flipping the backend and assert assign() restores it. + """ + import holoviews as hv + + from uxarray.plot.utils import HoloviewsBackend + + original = matplotlib.get_backend() + try: + matplotlib.use("svg") # the "user's" backend before plotting + + # Simulate hv.extension("matplotlib") flipping the active backend to agg. + monkeypatch.setattr(hv.Store, "current_backend", "bokeh", raising=False) + monkeypatch.setattr(hv, "extension", lambda *a, **k: matplotlib.use("agg")) + + be = HoloviewsBackend() + be.assign("matplotlib") + + # assign() must restore the backend that was active before the switch. + assert matplotlib.get_backend() == "svg" + finally: + matplotlib.use(original) + + +def test_inline_backend_reactivated_via_shell(monkeypatch): + """Inside IPython, restoring an inline backend must re-run the shell's own + backend activation (the public equivalent of ``%matplotlib inline``), which + rebuilds the full display integration that ``hv.extension`` tore down. + + Simply calling ``mpl.use`` is not enough: it restores the backend name but + leaves last-line figure auto-display broken (see #1537 / #1538). + """ + import sys + import types + + from uxarray.plot.utils import HoloviewsBackend + + inline_backend = "module://matplotlib_inline.backend_inline" + calls = [] + + class FakeShell: + def enable_matplotlib(self, gui): + calls.append(("enable_matplotlib", gui)) + + ipython = types.ModuleType("IPython") + ipython.get_ipython = lambda: FakeShell() + + monkeypatch.setitem(sys.modules, "IPython", ipython) + monkeypatch.setattr( + matplotlib, "use", lambda backend: calls.append(("use", backend)) + ) + + be = HoloviewsBackend() + be.matplotlib_backend = inline_backend + be.reset_mpl_backend() + + # The module:// inline backend must be mapped to the "inline" gui name and + # restored through the shell, without falling back to mpl.use. + assert calls == [("enable_matplotlib", "inline")] + + +def test_reset_backend_falls_back_to_mpl_use_outside_ipython(monkeypatch): + """Outside IPython (get_ipython() is None), restoration falls back to + ``mpl.use`` with the stored backend.""" + import sys + import types + + from uxarray.plot.utils import HoloviewsBackend + + calls = [] + ipython = types.ModuleType("IPython") + ipython.get_ipython = lambda: None + monkeypatch.setitem(sys.modules, "IPython", ipython) + monkeypatch.setattr(matplotlib, "use", lambda backend: calls.append(("use", backend))) + + be = HoloviewsBackend() + be.matplotlib_backend = "svg" + be.reset_mpl_backend() + + assert calls == [("use", "svg")] + + def test_plot_with_features(gridpath, datasetpath): """ensure can render a multiplot layout with geoviews features. Regression test for issue #1542. diff --git a/test/test_remap_yac.py b/test/test_remap_yac.py index 51f2b9e33..3efcd7c11 100644 --- a/test/test_remap_yac.py +++ b/test/test_remap_yac.py @@ -140,6 +140,23 @@ def test_yac_bilinear_face_remap(gridpath): assert out.size == dest.n_face +def test_yac_dnn_face_remap(gridpath): + # distance-nearest-neighbour (YAC >= 3.15); the default CELL_AREA search + # distance requires a face-centered target. + mesh_path = gridpath("mpas", "QU", "mesh.QU.1920km.151026.nc") + uxds = ux.open_dataset(mesh_path, mesh_path) + dest = ux.open_grid(mesh_path) + + out = uxds["latCell"].remap( + destination_grid=dest, + remap_to="faces", + backend="yac", + yac_method="dnn", + ) + + assert out.size == dest.n_face + + def test_yac_bilinear_rejects_non_average_method(gridpath): mesh_path = gridpath("mpas", "QU", "mesh.QU.1920km.151026.nc") uxds = ux.open_dataset(mesh_path, mesh_path) diff --git a/uxarray/grid/arcs.py b/uxarray/grid/arcs.py index c6d367024..1d0c5592d 100644 --- a/uxarray/grid/arcs.py +++ b/uxarray/grid/arcs.py @@ -8,7 +8,7 @@ _normalize_xyz_scalar, ) from uxarray.grid.utils import _angle_of_2_vectors -from uxarray.utils.computing import diff_of_products, two_sum +from uxarray.utils.computing import accucross, two_sum # Tolerance used to classify orient3d results as zero. For double-precision # unit-vector inputs this covers rounding error in the compensated cross product. @@ -406,6 +406,24 @@ def _orient3d_on_sphere_value(a, b, q): ) +@njit(cache=True, inline="always") +def _normal_dot_value(nx_hi, ny_hi, nz_hi, nx_lo, ny_lo, nz_lo, q0, q1, q2): + """Compensated dot of an already-computed ``(hi, lo)`` normal with ``q``. + + This is the second half of :func:`_orient3d_on_sphere_value_xyz`, split out + so that callers holding a normal that is invariant across many queries (e.g. + the ray plane ``q x R`` in the spherical point-in-polygon kernel, which is + the same for every edge of a face) can compute the cross product once and + reuse it, paying only this dot per query. + """ + p0 = (nx_hi + nx_lo) * q0 + p1 = (ny_hi + ny_lo) * q1 + p2 = (nz_hi + nz_lo) * q2 + s, e = two_sum(p0, p1) + s, e2 = two_sum(s, p2) + return s + (e + e2) + + @njit(cache=True, inline="always") def _orient3d_on_sphere_value_xyz(a0, a1, a2, b0, b1, b2, q0, q1, q2): """Scalar-argument form of :func:`_orient3d_on_sphere_value`. @@ -413,15 +431,8 @@ def _orient3d_on_sphere_value_xyz(a0, a1, a2, b0, b1, b2, q0, q1, q2): Takes the nine vector components directly so hot loops can call it without materializing ``(3,)`` arrays. """ - x_hi, x_lo = diff_of_products(a1, b2, a2, b1) - y_hi, y_lo = diff_of_products(a2, b0, a0, b2) - z_hi, z_lo = diff_of_products(a0, b1, a1, b0) - p0 = (x_hi + x_lo) * q0 - p1 = (y_hi + y_lo) * q1 - p2 = (z_hi + z_lo) * q2 - s, e = two_sum(p0, p1) - s, e2 = two_sum(s, p2) - return s + (e + e2) + nx_hi, ny_hi, nz_hi, nx_lo, ny_lo, nz_lo = accucross(a0, a1, a2, b0, b1, b2) + return _normal_dot_value(nx_hi, ny_hi, nz_hi, nx_lo, ny_lo, nz_lo, q0, q1, q2) @njit(cache=True) diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py index f4f4d96ff..f79c6586a 100644 --- a/uxarray/grid/grid.py +++ b/uxarray/grid/grid.py @@ -1951,8 +1951,21 @@ def calculate_total_face_area( ------- Sum of area of all the faces in the mesh : float """ + # Default parameters match the cached ``face_areas`` property, which also + # preserves the equal-area values used for HEALPix grids; reuse it to avoid + # recomputation. Any non-default quadrature settings require a fresh + # geometric computation so the requested rule/order actually take effect. + if ( + quadrature_rule == "triangular" + and order == 4 + and not latitude_adjusted_area + ): + return np.sum(self.face_areas.values) - return np.sum(self.face_areas.values) + face_areas, _ = self._compute_face_areas( + quadrature_rule, order, latitude_adjusted_area + ) + return np.sum(face_areas) def compute_face_areas( self, diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py index eaebb827b..94f8974ec 100644 --- a/uxarray/grid/intersections.py +++ b/uxarray/grid/intersections.py @@ -8,8 +8,7 @@ from uxarray.utils.computing import ( _cdp2, _cdp4, - _sum_sq_c2, - _sum_sq_c3, + _sum_of_squares_c, acc_sqrt_re, accucross, accucross_pair, @@ -471,9 +470,9 @@ def _accux_constlat_scalar(a0, a1, a2, b0, b1, b2, const_z): Invalid inputs propagate as non-finite coordinates. """ nx_hi, ny_hi, nz_hi, nx_lo, ny_lo, nz_lo = accucross(a0, a1, a2, b0, b1, b2) - s2_hi, s2_lo = _sum_sq_c2(nx_hi, nx_lo, ny_hi, ny_lo) + s2_hi, s2_lo = _sum_of_squares_c((nx_hi, ny_hi), (nx_lo, ny_lo)) denom = s2_hi + s2_lo - s3_hi, s3_lo = _sum_sq_c3(nx_hi, nx_lo, ny_hi, ny_lo, nz_hi, nz_lo) + s3_hi, s3_lo = _sum_of_squares_c((nx_hi, ny_hi, nz_hi), (nx_lo, ny_lo, nz_lo)) zsq_hi, zsq_lo = two_prod(const_z, const_z) d_hi, d_lo = _cdp4(s3_hi, zsq_hi, s3_hi, zsq_lo, s3_lo, zsq_hi, s3_lo, zsq_lo) e_hi, e_lo = two_sum(s2_hi, -d_hi) diff --git a/uxarray/grid/point_in_face.py b/uxarray/grid/point_in_face.py index a1f0b988f..3199c97d2 100644 --- a/uxarray/grid/point_in_face.py +++ b/uxarray/grid/point_in_face.py @@ -7,8 +7,13 @@ from numba import njit, prange from uxarray.constants import INT_DTYPE, INT_FILL_VALUE -from uxarray.grid.arcs import on_minor_arc, orient3d_on_sphere +from uxarray.grid.arcs import ( + _normal_dot_value, + _PREDICATE_ZERO_TOL, + on_minor_arc, +) from uxarray.grid.utils import _get_cartesian_face_edge_nodes +from uxarray.utils.computing import accucross if TYPE_CHECKING: from numpy.typing import ArrayLike @@ -69,8 +74,27 @@ def _ray_endpoint(q): return r -@njit(cache=True) -def _counts_as_crossing(A, B, q, R): +@njit(cache=True, inline="always") +def _sign_from_value(v): + """ + Sign of a compensated orient3d value under the standard zero tolerance. + """ + if v > _PREDICATE_ZERO_TOL: + return _SIGN_POS + if v < -_PREDICATE_ZERO_TOL: + return _SIGN_NEG + return _SIGN_ZERO + + +@njit(cache=True, inline="always") +def _counts_as_crossing( + a0, a1, a2, + b0, b1, b2, + q0, q1, q2, + r0, r1, r2, + qr_x_hi, qr_y_hi, qr_z_hi, + qr_x_lo, qr_y_lo, qr_z_lo, +): """Return 1 if edge AB crosses the minor arc q->R, 0 if not, -1 if degenerate. An edge AB crosses ray q->R iff q and R lie on opposite sides of the great @@ -79,12 +103,19 @@ def _counts_as_crossing(A, B, q, R): side-of-plane tests. Returns -1 when R lies exactly on plane(AB), which signals the caller to perturb R and retry. """ - s_AB_q = orient3d_on_sphere(A, B, q) - s_AB_R = orient3d_on_sphere(A, B, R) + # A x B: computed once, reused by both the q-side and R-side tests below. + nx_hi, ny_hi, nz_hi, nx_lo, ny_lo, nz_lo = accucross(a0, a1, a2, b0, b1, b2) + s_AB_q = _sign_from_value( + _normal_dot_value(nx_hi, ny_hi, nz_hi, nx_lo, ny_lo, nz_lo, q0, q1, q2) + ) # q on great circle AB: already caught by edge-membership check; not a crossing. if s_AB_q == _SIGN_ZERO: return 0 + + s_AB_R = _sign_from_value( + _normal_dot_value(nx_hi, ny_hi, nz_hi, nx_lo, ny_lo, nz_lo, r0, r1, r2) + ) # R on great circle AB: degenerate ray, caller must perturb R. if s_AB_R == _SIGN_ZERO: return -1 @@ -95,20 +126,27 @@ def _counts_as_crossing(A, B, q, R): # q and R are strictly on opposite sides of plane(AB). # Now check whether the intersection of the two great circles falls # inside the minor arc A->B, i.e. A and B are on opposite sides of plane(qR). - s_qR_A = orient3d_on_sphere(q, R, A) - s_qR_B = orient3d_on_sphere(q, R, B) + s_qR_A = _sign_from_value( + _normal_dot_value(qr_x_hi, qr_y_hi, qr_z_hi, qr_x_lo, qr_y_lo, qr_z_lo, a0, a1, a2) + ) + s_qR_B = _sign_from_value( + _normal_dot_value(qr_x_hi, qr_y_hi, qr_z_hi, qr_x_lo, qr_y_lo, qr_z_lo, b0, b1, b2) + ) - # A or B on great circle qR: vertex lies exactly on the ray plane. - # Apply the half-edge rule: count the edge only if the other endpoint is - # strictly on the negative side, so adjacent edges sharing this vertex - # are not double-counted. - if s_qR_A == _SIGN_ZERO or s_qR_B == _SIGN_ZERO: - if s_qR_A == _SIGN_ZERO and s_qR_B == _SIGN_ZERO: - return 0 # entire edge coplanar with ray: degenerate - s_other = s_qR_B if s_qR_A == _SIGN_ZERO else s_qR_A - return 1 if s_other == _SIGN_NEG else 0 + # Common case: neither endpoint lies on the ray plane, so the edge counts + # iff A and B straddle it. + s_qR_prod = s_qR_A * s_qR_B + if s_qR_prod != 0: + return 1 if s_qR_prod < 0 else 0 - return 1 if s_qR_A != s_qR_B else 0 + # An endpoint lies exactly on the ray plane. Apply the half-edge rule there: + # count the edge only if the other endpoint is strictly on the negative + # side, so that the two edges meeting at such a vertex are not both counted. + if s_qR_A == _SIGN_ZERO: + if s_qR_B == _SIGN_ZERO: + return 0 # whole edge coplanar with the ray plane: degenerate + return 1 if s_qR_B == _SIGN_NEG else 0 + return 1 if s_qR_A == _SIGN_NEG else 0 @njit(cache=True) @@ -163,13 +201,26 @@ def _point_in_polygon_sphere(q, polygon): # When R hits a degenerate edge, nudge and restart from i=0 so that all # edges are counted with the same ray — a mid-loop nudge corrupts parity. R = _ray_endpoint(q) + q0, q1, q2 = q[0], q[1], q[2] for _retry in range(4): + # The ray plane q x R is the same for every edge of the face, so it is + # computed once per ray pass rather than twice per edge. + qr_x_hi, qr_y_hi, qr_z_hi, qr_x_lo, qr_y_lo, qr_z_lo = accucross( + q0, q1, q2, R[0], R[1], R[2] + ) inside = False need_retry = False for i in range(n): A = polygon[i] B = polygon[(i + 1) % n] - c = _counts_as_crossing(A, B, q, R) + c = _counts_as_crossing( + A[0], A[1], A[2], + B[0], B[1], B[2], + q0, q1, q2, + R[0], R[1], R[2], + qr_x_hi, qr_y_hi, qr_z_hi, + qr_x_lo, qr_y_lo, qr_z_lo, + ) if c < 0: R[0] += 1e-7 R[1] -= 1e-7 diff --git a/uxarray/io/_scrip.py b/uxarray/io/_scrip.py index 3cc5b4d5c..ca4c513fb 100644 --- a/uxarray/io/_scrip.py +++ b/uxarray/io/_scrip.py @@ -9,6 +9,33 @@ from uxarray.grid.connectivity import _replace_fill_values +def _values_in_degrees(data_array): + """Return data_array.values, converting to degrees if necessary, + i.e., if data_array.attrs["units"] in ["radians", "radian", "rad"]. + + Parameters + ---------- + data_array : xr.DataArray + Array containing values to be extracted and converted to degrees. + Values assumed to be in degrees already unless data_array.attrs["units"] + is present and is one of ["radians", "radian", "rad"]. + + Returns + ------- + numpy.ndarray + data_array.values, converted to degrees if necessary. + """ + values = data_array.values + units = data_array.attrs.get("units", "").lower() + + # Check if units indicate radians + if units in ["radians", "radian", "rad"]: + return np.rad2deg(values) + + # Otherwise assume degrees (or no units means degrees for SCRIP) + return values + + def _to_ugrid(in_ds, out_ds): """If input dataset (``in_ds``) file is an unstructured SCRIP file, function will reassign SCRIP variables to UGRID conventions in output file @@ -21,8 +48,9 @@ def _to_ugrid(in_ds, out_ds): if any(key in in_ds for key in ["grid_imask", "grid_rank", "grid_area"]): # Create node_lon & node_lat variables from grid_corner_lat/lon # Turn latitude and longitude scrip arrays into 1D - corner_lat = in_ds["grid_corner_lat"].values.ravel() - corner_lon = in_ds["grid_corner_lon"].values.ravel() + # Convert to degrees if needed + corner_lat = _values_in_degrees(in_ds["grid_corner_lat"]).ravel() + corner_lon = _values_in_degrees(in_ds["grid_corner_lon"]).ravel() # Use Polars to find unique coordinate pairs df = pl.DataFrame({"lon": corner_lon, "lat": corner_lat}).with_row_count( @@ -62,14 +90,15 @@ def _to_ugrid(in_ds, out_ds): ) # Create face_lon & face_lat from grid_center_lat/lon + # Convert to degrees if needed out_ds[ugrid.FACE_COORDINATES[0]] = xr.DataArray( - in_ds["grid_center_lon"].values, + _values_in_degrees(in_ds["grid_center_lon"]), dims=[ugrid.FACE_DIM], attrs=ugrid.FACE_LON_ATTRS, ) out_ds[ugrid.FACE_COORDINATES[1]] = xr.DataArray( - in_ds["grid_center_lat"].values, + _values_in_degrees(in_ds["grid_center_lat"]), dims=[ugrid.FACE_DIM], attrs=ugrid.FACE_LAT_ATTRS, ) @@ -465,12 +494,17 @@ def _extract_single_grid( grid_corner_lat = grid_corner_lat.rename({corner_dim: "grid_corners"}) grid_corner_lon = grid_corner_lon.rename({corner_dim: "grid_corners"}) - grid_corner_lat = grid_corner_lat.copy() - grid_corner_lon = grid_corner_lon.copy() + # Convert to degrees if needed and create copies with correct units + grid_corner_lat_values = _values_in_degrees(grid_corner_lat) + grid_corner_lon_values = _values_in_degrees(grid_corner_lon) result = xr.Dataset() - result["grid_corner_lat"] = grid_corner_lat - result["grid_corner_lon"] = grid_corner_lon + result["grid_corner_lat"] = xr.DataArray( + grid_corner_lat_values, dims=grid_corner_lat.dims, attrs={"units": "degrees"} + ) + result["grid_corner_lon"] = xr.DataArray( + grid_corner_lon_values, dims=grid_corner_lon.dims, attrs={"units": "degrees"} + ) n_cells = grid_corner_lat.sizes["grid_size"] @@ -481,26 +515,40 @@ def _extract_single_grid( computed_lat_lon = None if center_lat and center_lat in ds: - center_lat_da = _stack_cell_dims( + stacked_lat = _stack_cell_dims( ds[center_lat], _resolve_cell_dims(metadata, ds[center_lat].dims), "grid_size", - ).copy() + ) + center_lat_da = xr.DataArray( + _values_in_degrees(stacked_lat), + dims=["grid_size"], + attrs={"units": "degrees"}, + ) else: if computed_lat_lon is None: computed_lat_lon = grid_center_lat_lon(result) - center_lat_da = xr.DataArray(computed_lat_lon[0], dims=["grid_size"]) + center_lat_da = xr.DataArray( + computed_lat_lon[0], dims=["grid_size"], attrs={"units": "degrees"} + ) if center_lon and center_lon in ds: - center_lon_da = _stack_cell_dims( + stacked_lon = _stack_cell_dims( ds[center_lon], _resolve_cell_dims(metadata, ds[center_lon].dims), "grid_size", - ).copy() + ) + center_lon_da = xr.DataArray( + _values_in_degrees(stacked_lon), + dims=["grid_size"], + attrs={"units": "degrees"}, + ) else: if computed_lat_lon is None: computed_lat_lon = grid_center_lat_lon(result) - center_lon_da = xr.DataArray(computed_lat_lon[1], dims=["grid_size"]) + center_lon_da = xr.DataArray( + computed_lat_lon[1], dims=["grid_size"], attrs={"units": "degrees"} + ) result["grid_center_lat"] = center_lat_da result["grid_center_lon"] = center_lon_da diff --git a/uxarray/plot/utils.py b/uxarray/plot/utils.py index e97919d6c..e5c8d3f8b 100644 --- a/uxarray/plot/utils.py +++ b/uxarray/plot/utils.py @@ -2,36 +2,68 @@ class HoloviewsBackend: - """Utility class to compare and set a HoloViews plotting backend for - visualization.""" + """Compare and set the HoloViews plotting backend.""" def __init__(self): self.matplotlib_backend = None def assign(self, backend: str): - """Assigns a backend for use with HoloViews visualization. - - Parameters - ---------- - backend : str - Plotting backend to use, one of 'matplotlib', 'bokeh' - """ - - if self.matplotlib_backend is None: - import matplotlib as mpl - - self.matplotlib_backend = mpl.get_backend() - + """Assign a HoloViews backend, one of 'matplotlib', 'bokeh'.""" if backend not in ["bokeh", "matplotlib", None]: raise ValueError( f"Unsupported backend. Expected one of ['bokeh', 'matplotlib'], but received {backend}" ) if backend is not None and backend != hv.Store.current_backend: - # only call hv.extension if it needs to be changed + import matplotlib as mpl + + # Capture the live backend now (not once at init) so a backend the + # user set later is what we restore. + self.matplotlib_backend = mpl.get_backend() hv.extension(backend) + if backend == "matplotlib": + # hv.extension("matplotlib") switches the active Matplotlib + # backend (e.g. to agg in Jupyter), breaking subsequent native + # matplotlib/xarray .plot() calls. HoloViews renders through + # hv.Store.current_backend, so restoring is safe. See #1537. + self.reset_mpl_backend() + def reset_mpl_backend(self): - """Resets the default backend for the ``matplotlib`` module.""" + """Restore the Matplotlib backend captured before the last switch. + + ``hv.extension("matplotlib")`` does not just change the active backend; + in an IPython/Jupyter kernel it also tears down the display integration + that auto-renders figures at the end of a cell. Simply calling + ``mpl.use`` puts the backend name back but leaves that integration + broken, so subsequent native ``matplotlib``/``xarray`` ``.plot()`` calls + silently produce no output unless ``plt.show()`` is called explicitly. + + Inside IPython we therefore re-run the shell's own backend activation + (the public equivalent of the ``%matplotlib`` magic), which rebuilds the + full integration in one step. Outside IPython we fall back to + ``mpl.use``. + """ + if self.matplotlib_backend is None: + return + + try: + from IPython import get_ipython + + shell = get_ipython() + except ImportError: + shell = None + + if shell is not None: + # Map the stored backend to the gui name enable_matplotlib expects. + gui = self.matplotlib_backend + if gui.startswith("module://") and "inline" in gui: + gui = "inline" + try: + shell.enable_matplotlib(gui) + return + except Exception: + pass + import matplotlib as mpl mpl.use(self.matplotlib_backend) diff --git a/uxarray/remap/accessor.py b/uxarray/remap/accessor.py index bb0f176ce..6415aa144 100644 --- a/uxarray/remap/accessor.py +++ b/uxarray/remap/accessor.py @@ -85,8 +85,10 @@ def nearest_neighbor( backend : {'uxarray', 'yac'}, default='uxarray' Remapping backend to use. When set to 'yac', requires YAC to be available on PYTHONPATH. - yac_method : {'nnn', 'average', 'conservative'}, optional + yac_method : {'nnn', 'dnn', 'average', 'conservative'}, optional YAC interpolation method. Defaults to 'nnn' when backend='yac'. + ``'dnn'`` (distance-nearest-neighbour, YAC >= 3.15) averages all + source points within a search distance of each target. yac_options : dict, optional YAC interpolation configuration options. @@ -206,10 +208,12 @@ def to_rectilinear( remapping before reshaping the result to latitude/longitude axes. The YAC backend uses YAC's rectilinear grid support directly and can be faster for large targets when YAC is installed. - yac_method : {'nnn', 'average', 'conservative'}, optional + yac_method : {'nnn', 'dnn', 'average', 'conservative'}, optional YAC interpolation method. When ``backend='yac'``, defaults to ``'nnn'`` because nearest-neighbor works for node-, edge-, and face-centered - source data. ``'conservative'`` requires face-centered source data. + source data. ``'dnn'`` (distance-nearest-neighbour) averages source + points within a search distance. ``'conservative'`` requires + face-centered source data. yac_options : dict, optional YAC interpolation configuration options forwarded to the selected YAC method. diff --git a/uxarray/remap/yac.py b/uxarray/remap/yac.py index 285a89a65..956ff175b 100644 --- a/uxarray/remap/yac.py +++ b/uxarray/remap/yac.py @@ -90,10 +90,11 @@ def _import_yac(): def _normalize_yac_method(yac_method: str | None) -> _YacOptions: if not yac_method: raise ValueError( - "backend='yac' requires yac_method to be set to 'nnn', 'average', or 'conservative'." + "backend='yac' requires yac_method to be set to 'nnn', 'dnn', " + "'average', or 'conservative'." ) method = yac_method.lower() - if method not in {"nnn", "average", "conservative"}: + if method not in {"nnn", "dnn", "average", "conservative"}: raise ValueError(f"Unsupported YAC method: {yac_method!r}") return _YacOptions(method=method, kwargs={}) @@ -232,6 +233,30 @@ def __init__( max_search_distance=yac_kwargs.get("max_search_distance", 0.0), scale=yac_kwargs.get("scale", 1.0), ) + elif yac_method == "dnn": + # Distance-nearest-neighbour (YAC >= 3.15): interpolate each target + # from all source points within a search distance around it. + weight_type = _coerce_enum( + yac_core.yac_interp_dnn_weight_type, + yac_kwargs.get("weight_type", yac_kwargs.get("reduction_type")), + ) + if weight_type is None: + weight_type = ( + yac_core.yac_interp_dnn_weight_type.YAC_INTERP_DNN_WEIGHT_DIST + ) + dnn_kwargs = {"weight_type": weight_type} + # search_distance_type default (CELL_AREA) requires a face-centered + # target; only forward it when explicitly provided. + if "search_distance_type" in yac_kwargs: + dnn_kwargs["search_distance_type"] = _coerce_enum( + yac_core.yac_interp_dnn_search_distance_type, + yac_kwargs["search_distance_type"], + ) + if "search_distance" in yac_kwargs: + dnn_kwargs["search_distance"] = yac_kwargs["search_distance"] + if "scale" in yac_kwargs: + dnn_kwargs["scale"] = yac_kwargs["scale"] + stack.add_dnn(**dnn_kwargs) elif yac_method == "average": reduction_type = _coerce_enum( yac_core.yac_interp_avg_weight_type, @@ -242,7 +267,7 @@ def __init__( yac_core.yac_interp_avg_weight_type.YAC_INTERP_AVG_ARITHMETIC ) stack.add_average( - reduction_type=reduction_type, + weight_type=reduction_type, partial_coverage=yac_kwargs.get("partial_coverage", False), ) elif yac_method == "conservative": diff --git a/uxarray/utils/computing.py b/uxarray/utils/computing.py index a0866bbd1..09ce972f1 100644 --- a/uxarray/utils/computing.py +++ b/uxarray/utils/computing.py @@ -81,45 +81,6 @@ def codegen(context, builder, signature, args): _FMA_INTRINSIC_OK = False -def _validate_fma() -> bool: - """Return True iff the FMA intrinsic compiles and is a bit-exact EFT.""" - if not _FMA_INTRINSIC_OK: - return False - try: - import numpy as _np - - @njit(cache=False) - def _tp_fma(a, b): - p = a * b - return p, _fma(a, b, -p) - - @njit(cache=False) - def _tp_vk(a, b): - p = a * b - f = 134217729.0 - a_hi = f * a - (f * a - a) - a_lo = a - a_hi - b_hi = f * b - (f * b - b) - b_lo = b - b_hi - e = a_lo * b_lo - (((p - a_hi * b_hi) - a_lo * b_hi) - a_hi * b_lo) - return p, e - - rng = _np.random.default_rng(20260101) - for _ in range(20000): - a = float(rng.standard_normal() * rng.integers(1, 1 << 20)) - b = float(rng.standard_normal() * rng.integers(1, 1 << 20)) - pf, ef = _tp_fma(a, b) - pv, ev = _tp_vk(a, b) - if pf != pv or (pf + ef) != (pv + ev): - return False - return True - except Exception: # pragma: no cover - return False - - -_HAS_FMA = _validate_fma() - - @njit(cache=True, inline="always") def two_sum(a, b): """Knuth's TwoSum: return (s, e) with s = fl(a + b) and s + e = a + b exactly. @@ -148,6 +109,14 @@ def two_sum(a, b): return s, e +if _FMA_INTRINSIC_OK: + + @njit(cache=True, inline="always") + def _two_prod_fma(a, b): + p = a * b + return p, _fma(a, b, -p) + + @njit(cache=True, inline="always") def _two_prod_veltkamp(a, b): """Portable TwoProd via Veltkamp splitting (no FMA dependency). @@ -166,13 +135,36 @@ def _two_prod_veltkamp(a, b): return p, e -if _HAS_FMA: +def _validate_fma(n_samples=2000) -> bool: + """Return True iff the FMA intrinsic compiles and is a bit-exact EFT.""" + if not _FMA_INTRINSIC_OK: + return False + try: + import numpy as _np - @njit(cache=True, inline="always") - def _two_prod_fma(a, b): - """TwoProd via a single hardware FMA: e = fma(a, b, -p).""" - p = a * b - return p, _fma(a, b, -p) + rng = _np.random.default_rng(20260101) + for _ in range(n_samples): + a = float(rng.standard_normal() * rng.integers(1, 1 << 20)) + b = float(rng.standard_normal() * rng.integers(1, 1 << 20)) + pf, ef = _two_prod_fma(a, b) + pv, ev = _two_prod_veltkamp(a, b) + # Compare the residuals *directly*. Do not compare ``pf + ef`` + # against ``pv + ev``: both sums round straight back to the product + # (|e| <= ulp(p)/2 by construction), so that predicate collapses to + # ``pf != pv`` -- a tautology, since both are fl(a*b). The exact + # residual is unique and representable, so a correct FMA and the + # Veltkamp split must agree bit-for-bit. + if pf != pv or ef != ev: + return False + return True + except Exception: # pragma: no cover + return False + + +_HAS_FMA = _validate_fma() + + +if _HAS_FMA: @njit(cache=True, inline="always") def two_prod(a, b): @@ -480,27 +472,52 @@ def _cdp4(a0, b0, a1, b1, a2, b2, a3, b3): @njit(cache=True, inline="always") -def _sum_sq_c2(h0, l0, h1, l1): - """Compensated sum of squares for 2 (hi, lo) pairs: h0²+l0²+h1²+l1². +def _fast_two_sum(a, b): + """Fast TwoSum: (x, y) with x = fl(a + b), x + y = a + b exactly. - Mirrors sum_of_squares_c from accusphgeom/numeric/eft.hpp, which - constructs lhs = rhs = [h0, l0, h1, l1] and calls compensated_dot_product. - Used to compute nx²+ny² accurately from the compensated normal (hi, lo). + Requires ``|a| >= |b|`` (or ``a == 0``) for the error term to be exact — + the callers here satisfy this because ``a`` is a running non-negative sum. """ - return _cdp4(h0, h0, l0, l0, h1, h1, l1, l1) + x = a + b + return x, (a - x) + b @njit(cache=True, inline="always") -def _sum_sq_c3(h0, l0, h1, l1, h2, l2): - """Compensated sum of squares for 3 (hi, lo) pairs: h0²+l0²+h1²+l1²+h2²+l2². +def _sum_non_neg(a_hi, a_lo, b_hi, b_lo): + """Add two non-negative compensated values (mirrors accusphgeom ``sum_non_neg``).""" + hh, h = two_sum(a_hi, b_hi) + d = h + (a_lo + b_lo) + return _fast_two_sum(hh, d) - Mirrors sum_of_squares_c from accusphgeom/numeric/eft.hpp, which - constructs lhs = rhs = [h0, l0, h1, l1, h2, l2] and calls a 6-term CDP. - We use _cdp8 with two zero-padding pairs (adding zero products). - Used to compute |n|² = nx²+ny²+nz² accurately from the compensated normal. + +@njit(cache=True, inline="always") +def _sum_of_squares_c(hi, lo): + """Compensated squared norm ``Σ (hi[i] + lo[i])²`` of a compensated vector. + + Faithful port of ``sum_of_squares_c`` from + accusphgeom/numeric/eft.hpp: a compensated ``Σ hi²`` plus the cross-term + correction ``2·Σ hi·li``. + + Parameters + ---------- + hi, lo : tuple of float + Equal-length tuples of the high and low parts of each vector component. + Numba specializes this per tuple length at compile time and keeps the + tuples register-resident (no allocation) — the direct analog of the C++ + template. Used for both nx²+ny² (denominator) and nx²+ny²+nz² (|n|²). """ - # lhs = rhs = [h0, l0, h1, l1, h2, l2, 0, 0] - return _cdp8(h0, l0, h1, l1, h2, l2, 0.0, 0.0, h0, l0, h1, l1, h2, l2, 0.0, 0.0) + n = len(hi) + s_hi = 0.0 + s_lo = 0.0 + for i in range(n): # compensated Σ hi² + ph, pl = two_prod(hi[i], hi[i]) + s_hi, s_lo = _sum_non_neg(s_hi, s_lo, ph, pl) + r_hi, r_lo = two_prod(hi[0], lo[0]) # accurate Σ hi·li + for i in range(1, n): + p, e = two_prod(hi[i], lo[i]) + r_hi, e2 = two_sum(r_hi, p) + r_lo += e + e2 + return _fast_two_sum(s_hi, (2.0 * (r_hi + r_lo)) + s_lo) @njit(cache=True, inline="always")