diff --git a/CHANGELOG.md b/CHANGELOG.md index 59a5c33..1630c43 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,31 @@ uses Semantic Versioning for public releases. ## Unreleased ### Fixed +- Nothing in an area result said how much of the sphere it was summed over. A + 5-degree mesh spanning 0-40E/0-40N with `sphere_radius: 6371000.0` returned + `total_area: 22936016715559.137 m^2` — 4.4967% of `4*pi*R^2` — with + `scientific_status: complete`, `physically_interpretable: true`, no warning + code, and a bare `postconditions: {status: not_evaluated, checks: []}`; the + identical call on a global mesh returned 1.0000 of the sphere with the same + status shape. `calculate_area` and `inspect_mesh` now both carry a + `mesh_coverage` block: `sphere_fraction`, `closed`, `euler_characteristic`, + `lon_extent`, `lat_extent`. A patch raises `MESH_NOT_GLOBAL` and drops to + `warning` but stays interpretable, because its total is a real physical + quantity and only the missing disclosure was wrong. The geometric and + topological halves are reported separately and are allowed to disagree: a + 1-degree structured global grid stops half a cell short of each pole, so it + reads `sphere_fraction: 0.999963` with `closed: false` and 720 boundary + edges, which is honest on both counts and not a regional patch. Counting + edge incidences is a Python loop — 1.43 s at 196,608 faces, 5.99 s at + 786,432 — so above 250,000 faces the topological half is skipped and + `closed` comes back `null` with `topology_skipped` giving the reason, + never `false`. +- An abstained postcondition now says why it abstained. The area identity + holds only on a closed mesh, so a regional result came back + `{status: not_evaluated, checks: []}` and the payload never distinguished + that from a deployment running `verdict_policy: off`. The block now carries + `not_evaluated_because` when the server can name a reason, read off + `mesh_coverage` so naming it costs no second traversal of the mesh. - The response contract described a payload the server does not send. It declared a top-level `physically_interpretable` boolean that no code path emits — every producer nests that verdict inside `scientific_status` — and diff --git a/docs/api.rst b/docs/api.rst index 920290a..08b2ec8 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -16,6 +16,10 @@ These modules contain the pure computation logic, separate from MCP and I/O. :members: :undoc-members: +.. automodule:: uxarray_mcp.domain.mesh_coverage + :members: + :undoc-members: + .. automodule:: uxarray_mcp.domain.variable :members: :undoc-members: diff --git a/docs/tools.md b/docs/tools.md index e591046..1db1f0c 100644 --- a/docs/tools.md +++ b/docs/tools.md @@ -312,6 +312,34 @@ Today `calculate_area` is the operation with a closed-form reference: the face areas of a closed mesh must sum to `4*pi*R^2`, or `4*pi` on a unit sphere. The check abstains — status `not_evaluated` — whenever it cannot be trusted: an open or regional mesh, a missing total, or an unreadable grid. +When the server can name the reason, the block carries it as +**`not_evaluated_because`**; a bare `not_evaluated` was indistinguishable +from a deployment running `verdict_policy: off`. + +## Mesh coverage + +`calculate_area` and `inspect_mesh` both report a **`mesh_coverage`** block: +`sphere_fraction`, `closed`, `euler_characteristic`, `lon_extent` and +`lat_extent`. A 5-degree mesh spanning 0-40E/0-40N returns `total_area: +22936016715559.137 m^2`, which is 4.4967% of `4*pi*R^2`; the same call on a +global mesh returns 1.0000 of the sphere, and before this block nothing in +either payload said which was which. A patch is warned about +(`MESH_NOT_GLOBAL`) rather than refused: its total is a real physical +quantity, and only the missing disclosure was wrong. + +The geometric and topological halves answer different questions and can +honestly disagree. A 1-degree structured global grid reads `sphere_fraction: +0.999963` with `closed: false` and `euler_characteristic: 0` — it stops half +a cell short of each pole, so it covers essentially the whole sphere and is +genuinely open. It is not warned about, because 3.7e-5 of the sphere is not +a regional patch, but the area identity still abstains and says why. + +Counting edge incidences is a Python loop over every face — 1.43 s on a +196,608-face HEALPix mesh, 5.99 s at the next zoom level — so above 250,000 +faces the topological half does not run. `closed` and +`euler_characteristic` come back `null` with `topology_skipped` giving the +reason, rather than a verdict nobody computed. `sphere_fraction` is +vectorized and is always reported. `calculate_area` also declares which sphere it measured on. UXarray computes face areas on the unit sphere and never applies `sphere_radius`, so a global diff --git a/src/uxarray_mcp/domain/__init__.py b/src/uxarray_mcp/domain/__init__.py index 2f9d5ed..f046374 100644 --- a/src/uxarray_mcp/domain/__init__.py +++ b/src/uxarray_mcp/domain/__init__.py @@ -6,6 +6,11 @@ from .area import compute_area_stats from .mesh import is_healpix_spec, load_dataset, load_grid, parse_healpix_zoom +from .mesh_coverage import ( + compute_mesh_coverage, + mesh_coverage_warning_codes, + mesh_is_closed, +) from .profile_coverage import ( compute_profile_coverage, profile_coverage_warning_codes, @@ -30,6 +35,9 @@ "is_healpix_spec", "parse_healpix_zoom", "compute_area_stats", + "compute_mesh_coverage", + "mesh_coverage_warning_codes", + "mesh_is_closed", "compute_profile_coverage", "profile_coverage_warning_codes", "compute_target_coverage", diff --git a/src/uxarray_mcp/domain/area.py b/src/uxarray_mcp/domain/area.py index cfba3c2..76e2f3e 100644 --- a/src/uxarray_mcp/domain/area.py +++ b/src/uxarray_mcp/domain/area.py @@ -21,6 +21,8 @@ from typing import Any +from .mesh_coverage import compute_mesh_coverage + #: UXarray's default when a grid declares nothing, and its unit-sphere basis. UNIT_SPHERE_RADIUS = 1.0 @@ -55,13 +57,18 @@ def compute_area_stats(grid: Any, sphere_radius: float | None = None) -> dict: if hasattr(face_areas, "attrs") and "units" in face_areas.attrs: area_units = face_areas.attrs["units"] + steradians = float(face_areas.sum()) stats = { - "total_area": float(face_areas.sum()), + "total_area": steradians, "mean_area": float(face_areas.mean()), "min_area": float(face_areas.min()), "max_area": float(face_areas.max()), "area_units": area_units, "n_face": int(grid.n_face), + # Attached before scaling, and measured on the unit sphere whatever + # radius is applied below: a total is only readable as global or + # regional next to the fraction of the sphere it was summed over. + "mesh_coverage": compute_mesh_coverage(grid, steradians=steradians), } radius, source = resolve_sphere_radius(grid, sphere_radius) return apply_sphere_radius(stats, radius, source) diff --git a/src/uxarray_mcp/domain/mesh_coverage.py b/src/uxarray_mcp/domain/mesh_coverage.py new file mode 100644 index 0000000..64666c5 --- /dev/null +++ b/src/uxarray_mcp/domain/mesh_coverage.py @@ -0,0 +1,242 @@ +"""How much of the sphere a mesh covers, and what shape it is. + +``calculate_area`` returns a regional patch's total the same way it returns +a global one. Measured on a 5-degree mesh spanning 0-40E/0-40N with +``sphere_radius=6371000.0``: ``total_area`` 22936016715559.137 m^2, which is +4.4967% of ``4*pi*R^2``, delivered with ``physically_interpretable: True``, +no warning codes, and ``postconditions: not_evaluated`` that never says why. +The same call on a global mesh returns 1.0000 of the sphere. Nothing in the +payload separated them. + +Two independent measurements are reported because they answer different +questions and can honestly disagree: + +``sphere_fraction`` + ``sum(face_areas) / (4*pi)`` on the unit sphere. Geometric: how much + surface is actually covered. Reported instead of the raw steradian sum + rather than alongside it -- the two differ by a mathematical constant, + and this block rides on every area and inspection result under a byte + budget (#83). +``closed`` / ``euler_characteristic`` + Topological: whether every edge is shared by exactly two faces, and + ``V - E + F``. + +A 1-degree structured global grid shows why both are needed. It reads +``sphere_fraction`` 0.999962 and ``closed`` False, ``euler_characteristic`` +0 -- ``Grid.from_structured`` stops its nodes at +/-89.5 there rather than +extending to the poles, so the mesh has two small polar holes. It covers +essentially the whole sphere and is genuinely open, and a single verdict +would have had to suppress one of those facts. The 2-degree grid of the +same family does reach the poles: 1.000000 and closed, ``euler`` 2. The +regional patch is a disk: ``euler`` 1. + +Cost is why the topology half is size-guarded. Counting edge incidences is +a Python loop over every face; on a 196,608-face HEALPix mesh it takes +1.43 s on top of 0.76 s for ``n_edge``, and the next zoom level -- 786,432 +faces -- costs 9.7 s for the pair. ``face_areas`` is vectorized and stays +under 0.05 s across all of these, so the geometric half is always +computed and the topological half abstains above the threshold, saying so +rather than reporting ``closed: false`` for a mesh nobody looked at. +""" + +from __future__ import annotations + +import math +from typing import Any + +import numpy as np + +#: How far ``sphere_fraction`` may sit from 1.0 and still count as global. +#: +#: Quadrature error is three orders of magnitude smaller than this -- a +#: 162-face global mesh integrates to 1.000002 -- so the slack is not for +#: numerical noise. It is for meshes that are global in every sense a +#: caller cares about but leave a pinhole somewhere: the 1-degree grid +#: above misses 3.8e-5 of the sphere at its poles. A mesh missing more +#: than 0.1% is missing something a caller would want named. +GLOBAL_COVERAGE_TOLERANCE = 1e-3 + +#: Face count above which the topology checks abstain rather than run. +#: +#: Set just above HEALPix zoom 7 (196,608 faces, 2.2 s for the pair) and +#: below zoom 8 (786,432 faces, 9.7 s). An inspection call that takes ten +#: seconds to report a boolean is not worth the boolean. +TOPOLOGY_MAX_FACES = 250_000 + +#: Decimal places used when matching node coordinates. Six is ~0.1 m on +#: Earth's surface, far below any mesh spacing we deal with, and coarse +#: enough to absorb the round-trip through NetCDF float64 text. +_COORD_DECIMALS = 6 + +#: Decimal places the reported numbers are rounded to. Full float64 repr +#: costs ~20 characters each on a block that is re-sent every turn, and +#: buys nothing a caller can use: 1e-6 of the sphere is 510 km^2, and 1e-6 +#: degree is ~0.1 m. +_REPORT_DECIMALS = 6 + + +def _canonical_node_ids(grid: Any) -> list[int]: + """Map nodes onto identity by position, not by index. + + A structured global grid stores the 0/360 seam twice and every pole + once per meridian, so counting edges on raw indices reports boundary + edges on a mesh that is geometrically closed. Merging nodes that sit + at the same point -- with all pole nodes collapsing to one, since + longitude is meaningless there -- makes the count reflect the surface + rather than the storage layout. + """ + lon = np.asarray(grid.node_lon, dtype=float) % 360.0 + lat = np.asarray(grid.node_lat, dtype=float) + seen: dict[str, int] = {} + ids: list[int] = [] + for x, y in zip(lon, lat): + if abs(abs(y) - 90.0) < 1e-9: + key = f"pole{y:+.1f}" + else: + key = ( + f"{round(x, _COORD_DECIMALS) % 360:.6f}_{round(y, _COORD_DECIMALS):.6f}" + ) + ids.append(seen.setdefault(key, len(seen))) + return ids + + +def mesh_is_closed(grid: Any) -> bool: + """True when every edge is shared by exactly two faces. + + A closed mesh is the precondition for the ``4*pi*R^2`` identity. The + cheap version of this test -- comparing ``n_edge`` against Euler's + formula -- is wrong on meshes with holes, so count edge incidences + directly. + """ + try: + connectivity = np.asarray(grid.face_node_connectivity) + node_ids = _canonical_node_ids(grid) + except Exception: # pragma: no cover - mocked grids in unit tests + return False + + n_node = len(node_ids) + incidence: dict[tuple[int, int], int] = {} + for face in connectivity: + nodes: list[int] = [] + for raw in face: + index = int(raw) + if not 0 <= index < n_node: + continue # fill value: a face with fewer nodes than the max + node = node_ids[index] + if not nodes or nodes[-1] != node: + nodes.append(node) + # A ring stored with a repeated first/last node is one edge, not two. + if len(nodes) > 1 and nodes[0] == nodes[-1]: + nodes.pop() + if len(nodes) < 3: + continue # degenerate after merging coincident nodes + for index, node in enumerate(nodes): + other = nodes[(index + 1) % len(nodes)] + key = (min(node, other), max(node, other)) + incidence[key] = incidence.get(key, 0) + 1 + if not incidence: + return False + return all(count == 2 for count in incidence.values()) + + +def compute_mesh_coverage( + grid: Any, + *, + steradians: float | None = None, +) -> dict[str, Any]: + """Measure how much of the sphere ``grid`` covers and what shape it is. + + Parameters + ---------- + grid : ux.Grid + Loaded UXarray grid. + steradians : float | None + ``sum(face_areas)`` on the unit sphere, when the caller has already + computed it. Passed in from ``compute_area_stats`` so the sum is not + paid for twice; recomputed here when absent. Not itself reported -- + it becomes ``sphere_fraction``. + + Returns + ------- + dict + ``sphere_fraction`` (geometric), ``closed`` and + ``euler_characteristic`` (topological, ``None`` when skipped), + ``lon_extent`` and ``lat_extent``. ``topology_skipped`` appears only + when the mesh was too large to check, carrying the reason, so a + ``None`` verdict is never mistaken for a negative one. + + ``lon_extent`` describes the mesh in the grid's own longitude + convention and is not a globality test: a global mesh stored on + [-180, 180] with 20-degree cells reads [-170, 170]. Use + ``sphere_fraction`` for that. + """ + coverage: dict[str, Any] = { + "sphere_fraction": None, + "closed": None, + "euler_characteristic": None, + "lon_extent": None, + "lat_extent": None, + } + + if steradians is None: + try: + steradians = float(np.asarray(grid.face_areas).sum()) + except Exception: # pragma: no cover - mocked grids in unit tests + steradians = None + if steradians is not None and math.isfinite(steradians): + coverage["sphere_fraction"] = round( + float(steradians) / (4.0 * math.pi), _REPORT_DECIMALS + ) + + try: + lon = np.asarray(grid.node_lon, dtype=float) + lat = np.asarray(grid.node_lat, dtype=float) + if lon.size and lat.size: + coverage["lon_extent"] = [ + round(float(lon.min()), _REPORT_DECIMALS), + round(float(lon.max()), _REPORT_DECIMALS), + ] + coverage["lat_extent"] = [ + round(float(lat.min()), _REPORT_DECIMALS), + round(float(lat.max()), _REPORT_DECIMALS), + ] + except Exception: # pragma: no cover - mocked grids in unit tests + pass + + try: + n_face = int(grid.n_face) + except Exception: # pragma: no cover - mocked grids in unit tests + return coverage + + if n_face > TOPOLOGY_MAX_FACES: + coverage["topology_skipped"] = ( + f"{n_face} faces exceeds the {TOPOLOGY_MAX_FACES}-face limit for " + "counting edge incidences; closure was not checked." + ) + return coverage + + try: + coverage["euler_characteristic"] = int(grid.n_node) - int(grid.n_edge) + n_face + except Exception: # pragma: no cover - mocked grids in unit tests + pass + coverage["closed"] = mesh_is_closed(grid) + return coverage + + +def mesh_coverage_warning_codes(coverage: dict[str, Any]) -> list[str]: + """Stable codes for a mesh that is not the whole sphere. + + Silent on a skipped topology check: not knowing whether a mesh is + closed is not evidence that it is open, and ``topology_skipped`` in + the block already says nobody looked. + """ + fraction = coverage.get("sphere_fraction") + if fraction is None: + return [] + if fraction < 1.0 - GLOBAL_COVERAGE_TOLERANCE: + return ["MESH_NOT_GLOBAL"] + if fraction > 1.0 + GLOBAL_COVERAGE_TOLERANCE: + # More surface than a sphere has means faces overlap or are stored + # twice. Quadrature cannot produce this at 1e-3. + return ["MESH_COVERAGE_EXCEEDS_SPHERE"] + return [] diff --git a/src/uxarray_mcp/postconditions.py b/src/uxarray_mcp/postconditions.py index 4a13305..7cf5da3 100644 --- a/src/uxarray_mcp/postconditions.py +++ b/src/uxarray_mcp/postconditions.py @@ -42,6 +42,8 @@ import os from typing import Any, Callable, Literal +from uxarray_mcp.domain.mesh_coverage import mesh_is_closed + #: Verdict policies, in order of decreasing generosity to the caller. VerdictPolicy = Literal["full", "reference_only", "off"] @@ -155,6 +157,7 @@ def evaluate_area_postconditions( The check abstains entirely when the mesh is not closed, because on an open regional mesh ``4*pi*R^2`` is not the right number and a failing verdict there would be the server being wrong, not the mesh. + ``area_identity_abstention`` turns that silence into a stated reason. """ total_area = result.get("total_area") if total_area is None or grid_loader is None: @@ -165,7 +168,12 @@ def evaluate_area_postconditions( except Exception: # pragma: no cover - a load failure is the caller's error return [] - if not mesh_is_closed(grid): + # ``mesh_coverage`` already counted edge incidences, and on a large mesh + # that is seconds rather than milliseconds. Fall back to counting again + # only for a remote worker on an older build, which sends no block. + coverage = result.get("mesh_coverage") or {} + closed = coverage.get("closed") if "closed" in coverage else mesh_is_closed(grid) + if not closed: return [] basis = result.get("area_basis") or {} @@ -213,96 +221,78 @@ def evaluate_area_postconditions( ] -#: Decimal places used when matching node coordinates. Six is ~0.1 m on -#: Earth's surface, far below any mesh spacing we deal with, and coarse -#: enough to absorb the round-trip through NetCDF float64 text. -_COORD_DECIMALS = 6 - +def area_identity_abstention(result: dict[str, Any]) -> str | None: + """Say why ``sum(face_areas) == 4*pi*R^2`` was not evaluated. -def _canonical_node_ids(grid: Any) -> list[int]: - """Map nodes onto identity by position, not by index. + Read off the ``mesh_coverage`` block rather than the grid, so naming + the reason costs no second traversal of a mesh that may have taken + seconds to traverse once. Returns ``None`` when the abstention has no + explanation this function can give -- an unloadable grid, or a result + from a worker old enough to send no coverage -- because inventing one + would be worse than the silence it replaces. - A structured global grid stores the 0/360 seam twice and every pole - once per meridian, so counting edges on raw indices reports boundary - edges on a mesh that is geometrically closed. Merging nodes that sit - at the same point -- with all pole nodes collapsing to one, since - longitude is meaningless there -- makes the count reflect the surface - rather than the storage layout. + Kept to one sentence on purpose: the block is re-sent on every later + turn of a conversation, so every word here is paid for repeatedly + (#83). """ - import numpy as np - - lon = np.asarray(grid.node_lon, dtype=float) % 360.0 - lat = np.asarray(grid.node_lat, dtype=float) - seen: dict[str, int] = {} - ids: list[int] = [] - for x, y in zip(lon, lat): - if abs(abs(y) - 90.0) < 1e-9: - key = f"pole{y:+.1f}" - else: - key = ( - f"{round(x, _COORD_DECIMALS) % 360:.6f}_{round(y, _COORD_DECIMALS):.6f}" - ) - ids.append(seen.setdefault(key, len(seen))) - return ids + coverage = result.get("mesh_coverage") + if not coverage: + return None + + if coverage.get("topology_skipped"): + return ( + "The 4*pi*R^2 identity holds only on a closed mesh and closure " + f"was not checked: {coverage['topology_skipped']}" + ) + if coverage.get("closed") is False: + fraction = coverage.get("sphere_fraction") + extent = ( + f" It covers {fraction:.4%} of the sphere." + if isinstance(fraction, (int, float)) + else "" + ) + return ( + "The 4*pi*R^2 identity holds only on a closed mesh, and this one " + f"has at least one boundary edge.{extent}" + ) -def mesh_is_closed(grid: Any) -> bool: - """True when every edge is shared by exactly two faces. - - A closed mesh is the precondition for the ``4*pi*R^2`` identity. The - cheap version of this test -- comparing ``n_edge`` against Euler's - formula -- is wrong on meshes with holes, so count edge incidences - directly. Meshes here are small enough for that to be free. - """ - try: - import numpy as np - - connectivity = np.asarray(grid.face_node_connectivity) - node_ids = _canonical_node_ids(grid) - except Exception: # pragma: no cover - mocked grids in unit tests - return False - - n_node = len(node_ids) - incidence: dict[tuple[int, int], int] = {} - for face in connectivity: - nodes: list[int] = [] - for raw in face: - index = int(raw) - if not 0 <= index < n_node: - continue # fill value: a face with fewer nodes than the max - node = node_ids[index] - if not nodes or nodes[-1] != node: - nodes.append(node) - # A ring stored with a repeated first/last node is one edge, not two. - if len(nodes) > 1 and nodes[0] == nodes[-1]: - nodes.pop() - if len(nodes) < 3: - continue # degenerate after merging coincident nodes - for index, node in enumerate(nodes): - other = nodes[(index + 1) % len(nodes)] - key = (min(node, other), max(node, other)) - incidence[key] = incidence.get(key, 0) + 1 - if not incidence: - return False - return all(count == 2 for count in incidence.values()) + return None def postcondition_block( checks: list[dict[str, Any]], policy: VerdictPolicy, + *, + not_evaluated_reason: str | None = None, ) -> dict[str, Any]: """Assemble the block attached to every analysis result. Present even when nothing was checked: #84's point is that an explicit ``not_evaluated`` costs almost nothing and stops a caller implying more confidence than the computation supports. + + ``not_evaluated`` on its own turned out to be half the disclosure. A + regional mesh came back with ``{"status": "not_evaluated", + "checks": []}`` because the area identity abstains on an open mesh, + and the payload never said that was why -- indistinguishable from a + deployment that had checking switched off. ``not_evaluated_because`` + carries the reason when the caller can be told one. """ if not checks or policy == "off": - return { + block: dict[str, Any] = { "status": STATUS_NOT_EVALUATED, "checks": [], "independent_verification": False, } + reason = ( + f"{VERDICT_POLICY_ENV.lower()}={policy}: no check was run." + if policy == "off" + else not_evaluated_reason + ) + if reason: + block["not_evaluated_because"] = reason + return block if policy == "reference_only": return { diff --git a/src/uxarray_mcp/response_contract.py b/src/uxarray_mcp/response_contract.py index f26e6b1..1f47940 100644 --- a/src/uxarray_mcp/response_contract.py +++ b/src/uxarray_mcp/response_contract.py @@ -117,6 +117,20 @@ def _field( ), ] +#: Shared by ``calculate_area`` and ``inspect_mesh``, which report the same +#: block. Declared once so the two cannot drift into describing it +#: differently. +_MESH_COVERAGE_DESCRIPTION = ( + "How much of the sphere the mesh spans and what shape it is: " + "`sphere_fraction` (geometric), `closed` and " + "`euler_characteristic` (topological, null when the mesh was too large " + "to check -- see `topology_skipped`), and `lon_extent`/`lat_extent` in " + "the grid's own convention. Face counts alone do not distinguish a " + "global mesh from a regional patch; `sphere_fraction` does. The two " + "halves can honestly disagree: a mesh with small polar holes reads " + "`sphere_fraction` 0.999962 and `closed` false." +) + #: Per-family declarations. Keyed by operation name so a caller can ask #: about exactly the call it is about to make. _CONTRACTS: dict[str, dict[str, Any]] = { @@ -142,6 +156,9 @@ def _field( "area.", required=False, ), + _field( + "mesh_coverage", "object", _MESH_COVERAGE_DESCRIPTION, required=False + ), ], }, "inspect_mesh": { @@ -150,6 +167,9 @@ def _field( _field("n_face", "integer", "Number of faces."), _field("n_node", "integer", "Number of nodes."), _field("n_edge", "integer", "Number of edges."), + _field( + "mesh_coverage", "object", _MESH_COVERAGE_DESCRIPTION, required=False + ), ], }, "calculate_zonal_mean": { diff --git a/src/uxarray_mcp/tools/frontdoor.py b/src/uxarray_mcp/tools/frontdoor.py index 21c44f9..e0a4774 100644 --- a/src/uxarray_mcp/tools/frontdoor.py +++ b/src/uxarray_mcp/tools/frontdoor.py @@ -11,10 +11,12 @@ from typing import Any from uxarray_mcp.domain.anomaly_coverage import anomaly_coverage_warning_codes +from uxarray_mcp.domain.mesh_coverage import mesh_coverage_warning_codes from uxarray_mcp.domain.profile_coverage import profile_coverage_warning_codes from uxarray_mcp.domain.subset_coverage import subset_coverage_warning_codes from uxarray_mcp.domain.temporal_coverage import temporal_coverage_warning_codes from uxarray_mcp.postconditions import ( + area_identity_abstention, evaluate_area_postconditions, postcondition_block, resolve_verdict_policy, @@ -367,6 +369,16 @@ def _finalize_analysis_result( status = "warning" warning_codes.append("AREA_UNITS_UNDECLARED") physically_interpretable = False + coverage_codes = mesh_coverage_warning_codes(result.get("mesh_coverage") or {}) + if coverage_codes: + # Warned, not refused, and interpretability is left alone: the + # area of a patch is a real physical quantity. 22936016715559 + # m^2 of Earth's surface between 0-40E and 0-40N is true, and + # calling it uninterpretable would be a second wrong answer. + # What was missing is that nothing said it was a patch, so a + # caller could read it as the planet. + status = "warning" + warning_codes.extend(coverage_codes) # Refuses by default when a declared precondition fails: raises # PreconditionRefusal unless the caller passed the override token. @@ -438,11 +450,16 @@ def _finalize_analysis_result( # check does run, whether the verdict comes with it is a policy. policy = resolve_verdict_policy(verdict_policy) post_checks: list[dict[str, Any]] = [] + abstention: str | None = None if operation == "calculate_area" and policy != "off": post_checks = evaluate_area_postconditions( result, _grid_loader(result), policy=policy ) - result["postconditions"] = postcondition_block(post_checks, policy) + if not post_checks: + abstention = area_identity_abstention(result) + result["postconditions"] = postcondition_block( + post_checks, policy, not_evaluated_reason=abstention + ) return result diff --git a/src/uxarray_mcp/tools/inspection.py b/src/uxarray_mcp/tools/inspection.py index 0195d91..44c5bd5 100644 --- a/src/uxarray_mcp/tools/inspection.py +++ b/src/uxarray_mcp/tools/inspection.py @@ -9,6 +9,7 @@ from uxarray_mcp.domain import ( compute_area_stats, + compute_mesh_coverage, compute_variable_info, compute_zonal_anomaly_stats, compute_zonal_mean_stats, @@ -39,6 +40,10 @@ def _inspect_mesh_local(file_path: str) -> Dict[str, Any]: - n_edge: Number of edges (boundaries between nodes) - n_max_face_nodes: Maximum number of nodes per face - file_size_mb: Size of the file in megabytes + - mesh_coverage: How much of the sphere the mesh spans and what + shape it is -- sphere_fraction, closed, euler_characteristic, + lon_extent, lat_extent. Counts alone do not distinguish a + global mesh from a regional patch. Example: >>> inspect_mesh("/path/to/mesh.nc") @@ -62,6 +67,7 @@ def _inspect_mesh_local(file_path: str) -> Dict[str, Any]: "n_edge": int(grid.n_edge), "n_max_face_nodes": int(grid.n_max_face_nodes), "file_size_mb": 0.0, + "mesh_coverage": compute_mesh_coverage(grid), }, tool="inspect_mesh", inputs={"file_path": file_path}, @@ -92,6 +98,7 @@ def _inspect_mesh_local(file_path: str) -> Dict[str, Any]: "n_edge": int(grid.n_edge), "n_max_face_nodes": int(grid.n_max_face_nodes), "file_size_mb": round(file_size_mb, 2), + "mesh_coverage": compute_mesh_coverage(grid), "recommended_next_steps": [ call("calculate_area", "grid_path"), call("plot_mesh", "grid_path"), diff --git a/src/uxarray_mcp/typed_results.py b/src/uxarray_mcp/typed_results.py index 6765e71..a0607fc 100644 --- a/src/uxarray_mcp/typed_results.py +++ b/src/uxarray_mcp/typed_results.py @@ -230,7 +230,8 @@ def declared_output_schemas() -> dict[str, dict[str, Any]]: "description": ( "Checks on the value after computing. A check may report " "'not_evaluated', which is an explicit abstention rather " - "than a pass." + "than a pass; when the server can say why it abstained, " + "'not_evaluated_because' carries the reason." ), "additionalProperties": True, }, diff --git a/tests/test_mesh_coverage.py b/tests/test_mesh_coverage.py new file mode 100644 index 0000000..28d4347 --- /dev/null +++ b/tests/test_mesh_coverage.py @@ -0,0 +1,267 @@ +"""A regional total must not look like a global one (#33). + +``calculate_area`` reported a patch the same way it reported the planet. +Measured on a 5-degree mesh spanning 0-40E/0-40N with +``sphere_radius=6371000.0`` passed as an argument: ``total_area`` +22936016715559.137 m^2, which is 4.4967% of ``4*pi*R^2``, delivered with +``scientific_status {'status': 'complete', 'physically_interpretable': +True, 'warning_codes': []}`` and ``postconditions {'status': +'not_evaluated', 'checks': []}``. The identical call on a global mesh +returned 1.0000 of the sphere with the same status. Nothing in either +payload separated them, and the abstention -- correct in itself, since +``4*pi*R^2`` does not hold on an open mesh -- never said why it abstained. + +The tests here compare the two payloads rather than checking one in +isolation, because "these two are indistinguishable" was the defect. +""" + +from __future__ import annotations + +import math +import warnings + +import numpy as np +import pytest +import uxarray as ux + +from uxarray_mcp.domain import mesh_coverage as coverage_module +from uxarray_mcp.domain.mesh_coverage import ( + compute_mesh_coverage, + mesh_coverage_warning_codes, +) +from uxarray_mcp.response_contract import describe_response_contract +from uxarray_mcp.tools.frontdoor import run_analysis + +#: The radius the measurement above was taken on. +EARTH_RADIUS_M = 6371000.0 + + +def _write(tmp_path, name, lon, lat): + grid = ux.Grid.from_structured(lon=lon, lat=lat) + path = tmp_path / f"{name}.nc" + grid.to_xarray().to_netcdf(path) + return str(path) + + +@pytest.fixture +def patch_grid_file(tmp_path): + """The 5-degree 0-40E/0-40N patch from the measurement above.""" + return _write(tmp_path, "patch", np.arange(0, 41, 5.0), np.arange(0, 41, 5.0)) + + +@pytest.fixture +def global_grid_file(tmp_path): + """A closed 20-degree global mesh, for the side-by-side comparison.""" + return _write(tmp_path, "global", np.arange(0, 360, 20.0), np.arange(-80, 81, 20.0)) + + +@pytest.fixture +def polar_hole_grid_file(tmp_path): + """Global in longitude, 1-degree in latitude, and open at both poles. + + ``Grid.from_structured`` extends its cells half a step past the outermost + coordinate, so latitudes running -89 to 89 in steps of 1 stop at +/-89.5 + and leave a small cap uncovered at each pole. The mesh covers 0.999963 of + the sphere and has 720 boundary edges: geometrically complete, + topologically open. A single "is this global" verdict would have had to + suppress one of those two facts. + """ + return _write(tmp_path, "polar", np.arange(0, 360, 20.0), np.arange(-89, 90, 1.0)) + + +def _area(grid_file, **kwargs): + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + return run_analysis( + operation="calculate_area", + grid_path=grid_file, + sphere_radius=EARTH_RADIUS_M, + **kwargs, + ) + + +def _inspect(grid_file): + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + return run_analysis(operation="inspect_mesh", grid_path=grid_file) + + +class TestAPatchIsDistinguishableFromThePlanet: + def test_the_two_results_disagree_about_the_sphere( + self, state_dir, patch_grid_file, global_grid_file + ): + """The measurement that motivated this, asserted directly.""" + patch = _area(patch_grid_file) + whole = _area(global_grid_file) + + assert patch["mesh_coverage"]["sphere_fraction"] == pytest.approx( + 0.044967, abs=1e-6 + ) + assert whole["mesh_coverage"]["sphere_fraction"] == pytest.approx(1.0, abs=1e-3) + + # The totals themselves stay what they always were: correct sums over + # the faces each mesh actually has. + assert patch["total_area"] == pytest.approx(22936016715559.137, rel=1e-9) + assert whole["total_area"] == pytest.approx( + 4 * math.pi * EARTH_RADIUS_M**2, rel=1e-4 + ) + + def test_only_the_patch_is_warned_about( + self, state_dir, patch_grid_file, global_grid_file + ): + patch = _area(patch_grid_file) + whole = _area(global_grid_file) + assert "MESH_NOT_GLOBAL" in patch["scientific_status"]["warning_codes"] + assert "MESH_NOT_GLOBAL" not in whole["scientific_status"]["warning_codes"] + + def test_the_patch_total_stays_interpretable(self, state_dir, patch_grid_file): + """Warned, not refused, and not demoted. + + 22936016715559 m^2 of Earth's surface between 0-40E and 0-40N is a + true physical quantity. Flipping ``physically_interpretable`` to + False would be a second wrong answer on top of the first: the number + was never the problem, the missing disclosure was. + """ + patch = _area(patch_grid_file) + assert patch["outcome"] == "complete" + assert patch["scientific_status"]["physically_interpretable"] is True + + def test_the_abstention_says_why( + self, state_dir, patch_grid_file, global_grid_file + ): + """``not_evaluated`` alone was half the disclosure.""" + patch = _area(patch_grid_file)["postconditions"] + assert patch["status"] == "not_evaluated" + reason = patch["not_evaluated_because"] + assert "closed" in reason + assert "4.4967%" in reason + + # A mesh the identity does hold on gets the check, not a reason. + whole = _area(global_grid_file)["postconditions"] + assert whole["status"] == "checked" + assert "not_evaluated_because" not in whole + + +class TestGeometryAndTopologyAreReportedSeparately: + def test_a_mesh_can_be_open_and_still_cover_the_sphere( + self, state_dir, polar_hole_grid_file + ): + block = _inspect(polar_hole_grid_file)["mesh_coverage"] + assert block["sphere_fraction"] == pytest.approx(0.999963, abs=1e-5) + assert block["closed"] is False + # Two disks removed from a sphere: 2 - 2 = 0. + assert block["euler_characteristic"] == 0 + + def test_the_euler_characteristic_separates_a_sphere_from_a_disk( + self, state_dir, patch_grid_file, global_grid_file + ): + assert _inspect(global_grid_file)["mesh_coverage"]["euler_characteristic"] == 2 + assert _inspect(patch_grid_file)["mesh_coverage"]["euler_characteristic"] == 1 + + def test_a_near_global_open_mesh_is_not_warned_about( + self, state_dir, polar_hole_grid_file + ): + """Missing 3.7e-5 of the sphere is not a regional patch. + + The warning exists so a caller does not read a patch total as the + planet. Firing it on a mesh with pinhole polar caps would spend the + caller's attention on a rounding difference. + """ + result = _area(polar_hole_grid_file) + assert "MESH_NOT_GLOBAL" not in result["scientific_status"]["warning_codes"] + + def test_the_identity_still_abstains_on_it_and_says_so( + self, state_dir, polar_hole_grid_file + ): + """Not warned is not the same as verified. + + ``4*pi*R^2`` genuinely does not hold on a mesh with holes, however + small, so the check abstains -- and now names the boundary rather + than leaving an empty ``not_evaluated``. + """ + block = _area(polar_hole_grid_file)["postconditions"] + assert block["status"] == "not_evaluated" + assert "boundary edge" in block["not_evaluated_because"] + + +class TestTopologyAbstainsRatherThanGuessing: + def test_a_large_mesh_skips_closure_and_says_it_did( + self, state_dir, monkeypatch, global_grid_file + ): + """``closed: null`` must never read as ``closed: false``. + + Counting edge incidences is a Python loop over every face: 1.43 s on + a 196,608-face HEALPix mesh and 5.99 s at the next zoom level, on top + of 0.76 s and 3.68 s for ``n_edge``. Above the limit the check does + not run, and the block says that instead of reporting a verdict + nobody computed. + """ + monkeypatch.setattr(coverage_module, "TOPOLOGY_MAX_FACES", 1) + block = _inspect(global_grid_file)["mesh_coverage"] + assert block["closed"] is None + assert block["euler_characteristic"] is None + assert "162 faces" in block["topology_skipped"] + # The geometric half is vectorized and stays. + assert block["sphere_fraction"] == pytest.approx(1.0, abs=1e-3) + + def test_a_skipped_check_raises_no_warning_code(self): + """Not knowing is not evidence of a defect.""" + assert mesh_coverage_warning_codes({"sphere_fraction": None}) == [] + assert mesh_coverage_warning_codes({}) == [] + + def test_the_abstention_reason_reaches_the_postcondition_block( + self, state_dir, monkeypatch, global_grid_file + ): + monkeypatch.setattr(coverage_module, "TOPOLOGY_MAX_FACES", 1) + block = _area(global_grid_file)["postconditions"] + assert block["status"] == "not_evaluated" + assert "closure was not checked" in block["not_evaluated_because"] + + def test_overlapping_faces_get_their_own_code(self): + """More surface than a sphere has is not a coverage shortfall. + + Quadrature error is ~2e-6 on a coarse global mesh, three orders + below the tolerance, so a fraction above 1.001 means faces are + stored twice or overlap -- a different defect from a patch. + """ + assert mesh_coverage_warning_codes({"sphere_fraction": 1.05}) == [ + "MESH_COVERAGE_EXCEEDS_SPHERE" + ] + + +class TestInspectMeshDescribesWithoutJudging: + def test_inspect_mesh_carries_the_block(self, state_dir, patch_grid_file): + """Counts alone never said which part of the sphere they counted.""" + result = _inspect(patch_grid_file) + assert result["mesh_coverage"]["sphere_fraction"] == pytest.approx( + 0.044967, abs=1e-6 + ) + assert result["mesh_coverage"]["lat_extent"] == [-2.5, 42.5] + + def test_inspecting_a_regional_mesh_is_not_a_warning( + self, state_dir, patch_grid_file + ): + """The code fires where a number is claimed, not where one is described. + + ``inspect_mesh`` reports what a mesh is. A regional mesh is not a + problem to be flagged; it becomes one only when a scalar total is + presented as if it covered the planet. + """ + status = _inspect(patch_grid_file)["scientific_status"] + assert status["warning_codes"] == [] + assert status["status"] == "complete" + + +class TestTheBlockIsDeclared: + @pytest.mark.parametrize("operation", ["calculate_area", "inspect_mesh"]) + def test_the_contract_declares_mesh_coverage(self, operation): + contract = describe_response_contract(operation) + declared = {field["name"] for field in contract["fields"]} + assert "mesh_coverage" in declared + assert "mesh_coverage" not in contract["required"] + + def test_a_mocked_grid_yields_nulls_rather_than_an_exception(self): + """A worker sending something grid-shaped must not crash the block.""" + block = compute_mesh_coverage(object()) + assert block["sphere_fraction"] is None + assert block["closed"] is None diff --git a/tests/test_payload_budget.py b/tests/test_payload_budget.py index b716f0c..e24ffb9 100644 --- a/tests/test_payload_budget.py +++ b/tests/test_payload_budget.py @@ -43,8 +43,17 @@ #: length of the temporary file paths echoed back in ``_provenance.inputs``. #: Lowered across the board when #83 removed the caller paths that #: ``recommended_next_steps`` used to interpolate into every suggestion. +#: Raised from 1150 and 1800 for the ``mesh_coverage`` block (#33): 151 +#: bytes on each. Before it, a 5-degree mesh spanning 0-40E/0-40N returned +#: `total_area: 22936016715559.137 m^2` -- 4.4967% of `4*pi*R^2` -- with +#: `physically_interpretable: True`, no warning code, and a bare +#: `postconditions: not_evaluated`, and the identical call on a global mesh +#: returned 1.0000 of the sphere. Nothing in either payload said which was +#: which. The block was trimmed to earn those bytes: the raw steradian sum +#: was dropped (it is `4*pi * sphere_fraction`) and the floats rounded to +#: 1e-6, which is 510 km^2 of sphere and 0.1 m of arc. RESULT_BYTE_BUDGETS = { - "inspect_mesh": 1150, + "inspect_mesh": 1300, # Kept above inspect_mesh for the postcondition block (#84/#90): ~440 # bytes that took correct verification answers from 11/20 to 20/20 in # the study, which is the one payload increase we have evidence for. @@ -52,7 +61,7 @@ # bytes. Before it, a global mesh returned `total_area: 12.566371` -- # 4*pi steradians -- with `area_units: null` and no warning code, and did # so even on a grid whose file declared `sphere_radius: 6371000.0`. - "calculate_area": 1800, + "calculate_area": 1950, "inspect_variable": 1700, # Raised from 2050 for the bin-coverage block and its precondition (#23), # most of it the repair text. Before it, a regional mesh asked for bands