diff --git a/COMPAT.md b/COMPAT.md index ae888321..6b1ef1e8 100644 --- a/COMPAT.md +++ b/COMPAT.md @@ -54,6 +54,13 @@ unchanged, but the new path is preferred going forward. | Fold `mask_fn`, custom fold providers (`get_deformed_orientation`), `directional_regularisation` direction callables and arrays, and the P2 `steepness_wtfunc` received or returned local (computational-frame) points and directions | They receive world points and give world directions; the interpolator converts to its local frame once. With an identity frame nothing changes (spec 005, stage S2) | 2026-11-04 | | `copy()` returned an interpolator with no frame; a transform change after constraints existed gave a stale result | New payloads carry `format_version: 2`, `constraint_space`, `transform` and `bounding_box`. `copy()` keeps the frame. A transform change after constraints exist raises `StaleFrameError` at the next set, evaluate or save (spec 005, stage S3). Loading a legacy payload with no choice is deprecated: see the deprecation table | 2026-11-05 | | `GeologicalFeature.get_data()` (used by `Loop2DView.add_data` and `ModelExporter.save`) returned locations moved by the local origin and rotated a second time; the intrusion classification evaluated the model at such points | Both use the stored world locations. With an identity frame nothing changes (spec 005, stage S4) | 2026-11-05 | +| `GeologicalFeature.evaluate_gradient_misfit()` changed the stored constraint arrays in place, compared a local vector with a world vector and read one column. Its value was a dot product with an unnormalised model gradient | The name and the public shape stay. The method changes no stored array. It gives the cosine between the world data direction and the model gradient for each row, so the value differs from the old value. Use `evaluate_gradient_cosines()` for a 1-D array. Two families of different length give one row of shape `(1, N1 + N2)`; before, this raised `ValueError` (spec 005, stage S5) | 2026-11-05 | +| `evaluate_gradient_misfit()` and `evaluate_value_misfit()` rebuilt a stale model | Both methods write one `UserWarning` and return NaN. They do not rebuild. Build the model first. The new public method `GeologicalFeature.verification_report()` gives the constraint diagnostics (residuals and exclusions) and does not build the model. A PLI model ignores gradient data, so its report shows an `unattributed` record for the gradient rows (spec 005, stage S5) | 2026-11-05 | +| `evaluate_value_misfit()` divided by the field range, with a division by zero for a zero range | It keeps the normalisation by the field range. A zero or not finite range gives the absolute difference. `evaluate_value_residuals()` gives the unnormalised 1-D result (spec 005, stage S5) | 2026-11-05 | +| Assignment of `nsteps` or `step_vector` on a 3D structured grid, and `set_nelements` | The extent stays the same. Assignment gives a new step, not a product rule. `nsteps` must be whole and at least 2. `set_nelements` keeps `maximum`. Grid steps and solved values can change slightly (spec 005, stage S5) | 2026-11-05 | +| Assignment of `origin`, `nsteps` or `step_vector` on a `RectilinearGrid` | It raises `NotImplementedError`. The node arrays are fixed (spec 005, stage S5) | 2026-11-05 | +| The 2D structured grid: `nsteps_cells`, `n_cell_x` and `n_cell_y` were separate values, and `maximum` could be assigned | The grid derives `nsteps_cells`, `n_cell_x` and `n_cell_y` from `nsteps`. `maximum` is read-only (assignment raises `AttributeError`) and keeps its legacy value (spec 005, stage S5) | 2026-11-05 | +| `ConstraintDiagnosticsReport` had no residual or exclusion fields | It has three new optional fields: `physical_residuals`, `assembled_residuals` and `exclusions`. With an identity frame and no resize, the value misfit is the same as before (spec 005, stage S5) | 2026-11-05 | The following stable methods remain part of the documented public API surface after the package extraction work and are covered by the snapshot-based contract test: add_onlap_unconformity, add_unconformity, create_and_add_domain_fault, create_and_add_fault, create_and_add_fold_frame, create_and_add_folded_fold_frame, create_and_add_folded_foliation, create_and_add_foliation, evaluate_model, evaluate_model_gradient, from_file, get_fault_surfaces, get_feature_by_name, get_stratigraphic_surfaces, rescale, save, scale. diff --git a/LoopStructural/modelling/features/_geological_feature.py b/LoopStructural/modelling/features/_geological_feature.py index 22a17bf0..5c438ccc 100644 --- a/LoopStructural/modelling/features/_geological_feature.py +++ b/LoopStructural/modelling/features/_geological_feature.py @@ -6,7 +6,18 @@ from __future__ import annotations +import dataclasses +import warnings + import numpy as np +from loop_interpolation import DiscreteInterpolator +from loop_interpolation._diagnostics import ( + ConstraintDiagnosticsReport, + ResidualSummary, + build_exclusion_records, + direction_cosines, + summarise_residuals, +) from LoopStructural.utils.maths import gradient_from_tetrahedron, regular_tetraherdron_for_points @@ -205,12 +216,12 @@ def evaluate_gradient( while not resolved: resolved = True for f in self.faults: - v = ( + side = ( f[0] .evaluate_value(tetrahedron.reshape(-1, 3), fillnan='nearest') .reshape(tetrahedron.shape[0], 4) ) - flag = np.logical_or(np.all(v > 0, axis=1), np.all(v < 0, axis=1)) + flag = np.logical_or(np.all(side > 0, axis=1), np.all(side < 0, axis=1)) if np.any(~flag): logger.warning( f"Points are too close to fault {f[0].name}. Refining the tetrahedron" @@ -237,55 +248,285 @@ def evaluate_gradient( logger.info(f'Gradient calculated for {self.name}') return v - def evaluate_gradient_misfit(self): + @staticmethod + def _first_stale(builder, label: str) -> str | None: + """Find the first stale builder or interpolator below a builder. + + The method only reads flags. It follows the rule of + ``BaseBuilder.up_to_date``. A structural frame builder or a fault + builder has no flag of its own, so the method reads its component + builders. It also reads the faults of each builder. + + Parameters + ---------- + builder : object or None + Builder to check. + label : str + Name of the builder in the text of the result. + + Returns + ------- + str or None + Text that names the first stale object, or None if none is stale. """ + if builder is None: + return None + for fault in getattr(builder, "faults", None) or []: + fault_builder = getattr(fault, "builder", fault) + fault_name = getattr(fault, "name", "a fault") + reason = GeologicalFeature._first_stale(fault_builder, f"the fault {fault_name}") + if reason is not None: + return reason + components = getattr(builder, "builders", None) + if components is not None: + for index, component in enumerate(components): + name = getattr(component, "name", index) + reason = GeologicalFeature._first_stale(component, f"{label} ({name})") + if reason is not None: + return reason + return None + if getattr(builder, "_up_to_date", True) is False: + return f"The builder of {label} is stale. Build the model first." + interpolator = getattr(builder, "_interpolator", None) + if getattr(interpolator, "up_to_date", True) is False: + return f"The interpolator of {label} is stale. Build the model first." + return None + + def _stale_reason(self) -> str | None: + """Find the first stale object of the model. + + The method only reads flags. It does not build the model. Returns ------- - misfit : np.array(N,dtype=double) - dot product between interpolated gradient and constraints + str or None + Text that names the first stale object and says to build the model + first, or None if no object is stale. """ - if self.builder is not None: - self.builder.up_to_date() - grad = self.interpolator.get_gradient_constraints() - norm = self.interpolator.get_norm_constraints() - - dot = [] - if grad.shape[0] > 0: - grad /= np.linalg.norm(grad, axis=1)[:, None] - positions = grad[:, :3] - if self.interpolator.bounding_box is not None: - positions = self.interpolator.bounding_box.reproject(positions) - model_grad = self.evaluate_gradient(positions) - dot.append(np.einsum("ij,ij->i", model_grad, grad[:, :3:6]).tolist()) - - if norm.shape[0] > 0: - norm /= np.linalg.norm(norm, axis=1)[:, None] - positions = norm[:, :3] - if self.interpolator.bounding_box is not None: - positions = self.interpolator.bounding_box.reproject(positions) - model_norm = self.evaluate_gradient(positions) - dot.append(np.einsum("ij,ij->i", model_norm, norm[:, :3:6])) - - return np.array(dot) - - def evaluate_value_misfit(self): + reason = self._first_stale(self.builder, self.name) + if reason is not None: + return reason + for fault in getattr(self, "faults", None) or []: + reason = self._first_stale( + getattr(fault, "builder", None), f"the fault {getattr(fault, 'name', 'a fault')}" + ) + if reason is not None: + return reason + if getattr(self.interpolator, "up_to_date", True) is False: + return f"The interpolator of {self.name} is stale. Build the model first." + return None + + def _constraint_row_counts(self) -> list: + """Count the stored rows of the gradient, normal and value constraints. + + Returns + ------- + list of int + Row counts of the gradient, normal and value arrays. """ + return [ + int(np.shape(self.interpolator.get_gradient_constraints())[0]), + int(np.shape(self.interpolator.get_norm_constraints())[0]), + int(np.shape(self.interpolator.get_value_constraints())[0]), + ] + + def _gradient_cosines_fresh(self) -> list: + """Compute the gradient and normal cosines without a stale check. + + The method copies the stored arrays. It does not change them. Returns ------- - misfit : np.array(N,dtype=double) - difference between interpolated scalar field and value constraints + list of numpy.ndarray + One array of cosines for each non-empty family (gradient first, + then normal). """ - if self.builder is not None: - self.builder.up_to_date() + bounding_box = self.interpolator.bounding_box + cosines = [] + for array in ( + self.interpolator.get_gradient_constraints(), + self.interpolator.get_norm_constraints(), + ): + if np.shape(array)[0] == 0: + continue + rows = np.array(array, dtype=float, copy=True) + positions = rows[:, :3] + data = rows[:, 3:6] + if bounding_box is not None: + positions = bounding_box.reproject(positions) + data = self.interpolator.local_to_world_vectors(data) + model = self.evaluate_gradient(positions) + cosines.append(direction_cosines(model, data)) + return cosines + + def _value_residuals_fresh(self) -> np.ndarray: + """Compute the value residuals without a stale check. - locations = self.interpolator.get_value_constraints() - positions = locations[:, :3] + Returns + ------- + numpy.ndarray + Absolute difference in field units, shape (N,). + """ + rows = np.array(self.interpolator.get_value_constraints(), dtype=float, copy=True) + positions = rows[:, :3] if self.interpolator.bounding_box is not None: positions = self.interpolator.bounding_box.reproject(positions) - diff = np.abs(locations[:, 3] - self.evaluate_value(positions)) - diff /= self.max() - self.min() + return np.abs(rows[:, 3] - self.evaluate_value(positions)) + + def _physical_residual_summaries(self) -> dict[str, ResidualSummary]: + """Summarise the physical residuals of the value, gradient and normal data. + + The method does not check if the model is stale. The caller must check. + + Returns + ------- + dict of str to ResidualSummary + One summary for each non-empty family. The value summary uses field + units. The gradient and normal summaries use degrees. + """ + summaries: dict[str, ResidualSummary] = {} + counts = self._constraint_row_counts() + if counts[2] > 0: + summaries["value"] = summarise_residuals( + "value", "physical", "field", self._value_residuals_fresh() + ) + # The cosine list has one array for each non-empty family, in this order. + names = [name for name, count in zip(("gradient", "normal"), counts[:2]) if count > 0] + for name, cosine in zip(names, self._gradient_cosines_fresh()): + angles = np.degrees(np.arccos(np.clip(cosine, -1, 1))) + summaries[name] = summarise_residuals(name, "physical", "degree", angles) + return summaries + + def verification_report(self) -> ConstraintDiagnosticsReport: + """Build a report of the constraints, residuals and exclusions. + + The method returns a new object. It does not store the report, build the + model or change the state of the model. + + Returns + ------- + ConstraintDiagnosticsReport + The report. The residual blocks are empty if the model is not + solved. On a stale model the method warns, leaves the residual + blocks empty and still gives the exclusions. + """ + base = self.interpolator._build_constraint_diagnostics_report() + physical: dict[str, ResidualSummary] = {} + assembled: dict[str, ResidualSummary] = {} + reason = self._stale_reason() + if reason is not None: + warnings.warn(reason, UserWarning, stacklevel=2) + elif getattr(self.interpolator, "c", None) is not None: + physical = self._physical_residual_summaries() + if hasattr(self.interpolator, "calculate_residual_for_constraints"): + residuals = self.interpolator.calculate_residual_for_constraints() + for name, residual in residuals.items(): + assembled[name] = summarise_residuals(name, "assembled", "assembled", residual) + excluded = getattr(self.builder, "excluded_rows", {}) if self.builder is not None else {} + rows_per_point = ( + {"gradient": 2} if isinstance(self.interpolator, DiscreteInterpolator) else None + ) + return dataclasses.replace( + base, + physical_residuals=physical, + assembled_residuals=assembled, + exclusions=build_exclusion_records(excluded, base.families, rows_per_point), + ) + + def evaluate_gradient_cosines(self) -> np.ndarray: + """Compute the cosine between the model gradient and the gradient data. + + The cosine uses world coordinates (see COMPAT.md). The gradient rows come + first and the normal rows come second. The method does not change the + model or the stored constraints. + + Returns + ------- + numpy.ndarray + Cosines, shape (N,). The array is empty if there are no rows. On a + stale model the method warns and returns NaN for each row. + """ + reason = self._stale_reason() + if reason is not None: + warnings.warn(reason, UserWarning, stacklevel=2) + n = sum(self._constraint_row_counts()[:2]) + return np.full(n, np.nan) + cosines = self._gradient_cosines_fresh() + if len(cosines) == 0: + return np.zeros(0) + return np.concatenate(cosines) + + def evaluate_value_residuals(self) -> np.ndarray: + """Compute the residual between the model field and the value data. + + The method does not change the model or the stored constraints. + + Returns + ------- + numpy.ndarray + Absolute difference in field units, shape (N,), in the order of the + value rows. On a stale model the method warns and returns NaN for + each row. + """ + reason = self._stale_reason() + if reason is not None: + warnings.warn(reason, UserWarning, stacklevel=2) + return np.full(self._constraint_row_counts()[2], np.nan) + return self._value_residuals_fresh() + + @staticmethod + def _misfit_shape(counts: list) -> tuple: + """Give the shape of the old gradient misfit array for the row counts.""" + if len(counts) == 0: + return (0,) + if len(counts) == 2 and counts[0] == counts[1]: + return (2, counts[0]) + return (1, sum(counts)) + + def evaluate_gradient_misfit(self) -> np.ndarray: + """Compute the cosine misfit of the gradient and normal constraints. + + The value is a cosine in world coordinates (see COMPAT.md). + + Returns + ------- + numpy.ndarray + Shape (1, N) for one family, (2, N) for two families of equal + length, (1, N1 + N2) for two families of different length and (0,) + with no rows. On a stale model the method warns and returns NaN in + the same shape. + """ + reason = self._stale_reason() + if reason is not None: + warnings.warn(reason, UserWarning, stacklevel=2) + counts = [n for n in self._constraint_row_counts()[:2] if n > 0] + return np.full(self._misfit_shape(counts), np.nan) + cosines = self._gradient_cosines_fresh() + if len(cosines) == 0: + return np.zeros(0) + if len(cosines) == 2 and len(cosines[0]) == len(cosines[1]): + return np.vstack(cosines) + return np.concatenate(cosines)[None, :] + + def evaluate_value_misfit(self) -> np.ndarray: + """Compute the value misfit. + + Returns + ------- + numpy.ndarray + Absolute difference between the model field and the value data, + shape (N,). The difference is divided by the range of the field + only if the range is finite and larger than zero. On a stale model + the method warns and returns NaN for each row. + """ + reason = self._stale_reason() + if reason is not None: + warnings.warn(reason, UserWarning, stacklevel=2) + return np.full(self._constraint_row_counts()[2], np.nan) + diff = self._value_residuals_fresh() + value_range = self.max() - self.min() + if np.isfinite(value_range) and value_range > 0: + diff = diff / value_range return diff def copy(self, name=None): diff --git a/ROADMAP.md b/ROADMAP.md index 1a24ccd3..335e1fc2 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -140,6 +140,18 @@ Removal is not before 2 minor releases after the first release that warns. Removal is also not before the known callers (the QGIS plugin) have moved. The library no longer calls either method. +### Misfit diagnostics of GeologicalFeature (spec 005, S5) +Decision date: 2026-11-05 (spec 005, S5, Q2). + +The misfit methods `evaluate_gradient_misfit` and `evaluate_value_misfit` +keep their names and their public shape. The gradient value is now a cosine +in world coordinates. + +A stale model gives a `UserWarning` and NaN. The model is not rebuilt. + +The new methods `evaluate_gradient_cosines` and `evaluate_value_residuals` +give the 1-D result. The notice is in `COMPAT.md`. + ## Stage sequence - [x] **Stage 0 — Planning infra.** This file, release/versioning/compat @@ -881,3 +893,4 @@ The library no longer calls either method. regressions. - **2026-11-05:** Spec 005 stage S3 implemented: saved interpolator frame, legacy-load rule, frame-keeping copies, stale-frame check. Gate result: root pytest 1978 passed, 5 skipped. - **2026-11-05:** Spec 005 stage S4 implemented: display, export and intrusion classification use stored world locations; rescale and scale deprecated. +- **2026-11-05:** Spec 005 stage S5 implemented: resize paths keep the extent, cell-centre guard tests, non-mutating world-correct misfit diagnostics, residual and exclusion report. diff --git a/packages/loop_common/src/loop_common/geometry/_structured_grid_2d.py b/packages/loop_common/src/loop_common/geometry/_structured_grid_2d.py index 7c48836b..16244673 100644 --- a/packages/loop_common/src/loop_common/geometry/_structured_grid_2d.py +++ b/packages/loop_common/src/loop_common/geometry/_structured_grid_2d.py @@ -1,6 +1,5 @@ """Pure 2D regular grid geometry: origin/nsteps/step_vector indexing.""" - import numpy as np from loop_common.logging import get_logger as getLogger @@ -20,14 +19,55 @@ def __init__(self, origin=None, nsteps=None, step_vector=None): nsteps = np.array([10, 10]) if step_vector is None: step_vector = np.ones(2) - self.nsteps = np.ceil(np.array(nsteps)).astype(int) + self.nsteps = nsteps self.step_vector = np.array(step_vector) self.origin = np.array(origin) - self.maximum = origin + self.nsteps * self.step_vector self.dim = 2 - self.nsteps_cells = self.nsteps - 1 - self.n_cell_x = self.nsteps[0] - 1 - self.n_cell_y = self.nsteps[1] - 1 + + @property + def nsteps(self) -> np.ndarray: + """Number of nodes on each axis. + + Returns + ------- + np.ndarray + Integer array of node counts. + """ + return self._nsteps + + @nsteps.setter + def nsteps(self, nsteps) -> None: + # Round up to whole nodes. All counts come from this value. + self._nsteps = np.ceil(np.asarray(nsteps)).astype(int) + + @property + def nsteps_cells(self) -> np.ndarray: + """Number of cells on each axis (`nsteps - 1`).""" + return self.nsteps - 1 + + @property + def n_cell_x(self): + """Number of cells on the x axis.""" + return self.nsteps[0] - 1 + + @property + def n_cell_y(self): + """Number of cells on the y axis.""" + return self.nsteps[1] - 1 + + @property + def maximum(self) -> np.ndarray: + """Legacy maximum corner (FR-028). + + Returns + ------- + np.ndarray + `origin + nsteps * step_vector`. + """ + # Legacy value (FR-028). It counts nodes as cells, so it is one + # step beyond the last node. `nodes` uses `nsteps_cells`. + # Do not change it. + return self.origin + self.nsteps * self.step_vector @property def nodes(self): @@ -88,7 +128,10 @@ def inside(self, pos: np.ndarray) -> np.ndarray: inside = np.ones(pos.shape[0]).astype(bool) for i in range(self.dim): inside *= pos[:, i] > self.origin[None, i] - inside *= pos[:, i] < self.origin[None, i] + self.step_vector[None, i] * self.nsteps_cells[None, i] + inside *= ( + pos[:, i] + < self.origin[None, i] + self.step_vector[None, i] * self.nsteps_cells[None, i] + ) return inside def check_position(self, pos: np.ndarray) -> np.ndarray: @@ -118,14 +161,20 @@ def neighbour_global_indexes(self, mask=None, **kwargs): if mask is None: mask = np.array([[-1, 0, 1, -1, 0, 1, -1, 0, 1], [1, 1, 1, 0, 0, 0, -1, -1, -1]]) neighbours = indexes[:, None, :] + mask[:, :, None] - return (neighbours[0, :, :] + self.nsteps[0, None, None] * neighbours[1, :, :]).astype(np.int64) + return (neighbours[0, :, :] + self.nsteps[0, None, None] * neighbours[1, :, :]).astype( + np.int64 + ) def cell_corner_indexes(self, cell_indexes: np.ndarray) -> np.ndarray: corner_indexes = np.zeros((cell_indexes.shape[0], 4, 2), dtype=np.int64) xcorner = np.array([0, 1, 0, 1]) ycorner = np.array([0, 0, 1, 1]) - corner_indexes[:, :, 0] = cell_indexes[:, None, 0] + corner_indexes[:, :, 0] + xcorner[None, :] - corner_indexes[:, :, 1] = cell_indexes[:, None, 1] + corner_indexes[:, :, 1] + ycorner[None, :] + corner_indexes[:, :, 0] = ( + cell_indexes[:, None, 0] + corner_indexes[:, :, 0] + xcorner[None, :] + ) + corner_indexes[:, :, 1] = ( + cell_indexes[:, None, 1] + corner_indexes[:, :, 1] + ycorner[None, :] + ) return corner_indexes def global_index_to_cell_index(self, global_index): diff --git a/packages/loop_common/src/loop_common/geometry/_structured_grid_3d.py b/packages/loop_common/src/loop_common/geometry/_structured_grid_3d.py index cf1c3514..30723952 100644 --- a/packages/loop_common/src/loop_common/geometry/_structured_grid_3d.py +++ b/packages/loop_common/src/loop_common/geometry/_structured_grid_3d.py @@ -1,4 +1,6 @@ -"""Pure 3D regular grid geometry: origin/nsteps/step_vector indexing.""" +"""Pure 3D regular grid geometry: origin/nsteps/step_vector indexing.""" + +from typing import Optional import numpy as np @@ -52,7 +54,8 @@ def set_nelements(self, nelements) -> int: ele_vol = box_vol / nelements step_vector = np.zeros(3) step_vector[:] = ele_vol ** (1.0 / 3.0) - nsteps = np.ceil((self.maximum - self.origin) / step_vector).astype(int) + # Keep the extent. Use at least 2 nodes on each axis. + nsteps = np.maximum(np.ceil((self.maximum - self.origin) / step_vector).astype(int), 2) self.nsteps = nsteps return self.n_elements @@ -70,9 +73,16 @@ def nsteps(self): @nsteps.setter def nsteps(self, nsteps): - change_factor = nsteps / self.nsteps - self._step_vector /= change_factor - self._nsteps = nsteps + nsteps = np.broadcast_to(np.asarray(nsteps), self._nsteps.shape) + if not np.all(np.isfinite(nsteps)) or np.any(nsteps != np.round(nsteps)): + raise ValueError("Cannot resize the grid. nsteps must be whole numbers.") + nsteps = nsteps.astype(int) + if np.any(nsteps < 2): + raise ValueError("Cannot resize the grid. nsteps must be at least 2.") + # Keep the extent: the new step is the extent divided by the new cell count. + extent = self.nsteps_cells * self._step_vector + self._step_vector = extent / (nsteps - 1) + self._nsteps = np.array(nsteps) self._on_geometry_change() @property @@ -97,13 +107,17 @@ def rotation_xy(self, rotation_xy): ) self._rotation_xy = validate_rigid_rotation(rotation_xy, 3) - def _cells_for_length(self, length: np.ndarray) -> np.ndarray: + def _cells_for_length( + self, length: np.ndarray, step_vector: Optional[np.ndarray] = None # noqa: FA100 + ) -> np.ndarray: """Return cell counts that cover the requested axis lengths. Parameters ---------- length : np.ndarray Requested support length along each grid axis. + step_vector : Optional[np.ndarray] + Step to use along each axis. The default is the current step. Returns ------- @@ -115,7 +129,9 @@ def _cells_for_length(self, length: np.ndarray) -> np.ndarray: ValueError If any requested cell count is negative. """ - cells = np.ceil(length / self.step_vector - 1e-9).astype(np.int64) + if step_vector is None: + step_vector = self.step_vector + cells = np.ceil(length / step_vector - 1e-9).astype(np.int64) cells[cells == 0] = 2 if np.any(cells < 0): logger.error( @@ -130,9 +146,11 @@ def step_vector(self): @step_vector.setter def step_vector(self, step_vector): - change_factor = step_vector / self._step_vector - newsteps = self._nsteps / change_factor - self._nsteps = np.ceil(newsteps).astype(int) + step_vector = np.array(step_vector, dtype=float) + # Keep the extent: the new cell count covers the old extent with the new step. + extent = self.nsteps_cells * self._step_vector + cells = self._cells_for_length(extent, step_vector=step_vector) + self._nsteps = cells + 1 self._step_vector = step_vector self._on_geometry_change() diff --git a/packages/loop_common/src/loop_common/supports/_3d_rectilinear_grid.py b/packages/loop_common/src/loop_common/supports/_3d_rectilinear_grid.py index 52e8abb1..1a251c45 100644 --- a/packages/loop_common/src/loop_common/supports/_3d_rectilinear_grid.py +++ b/packages/loop_common/src/loop_common/supports/_3d_rectilinear_grid.py @@ -16,6 +16,10 @@ logger = getLogger(__name__) +_RESIZE_MESSAGE = ( + "A rectilinear grid has fixed node arrays and cannot be resized. Build a new RectilinearGrid." +) + class RectilinearGrid(StructuredGrid): """A rectilinear (non-uniformly spaced) 3-D structured grid. @@ -78,6 +82,91 @@ def ynodes(self) -> np.ndarray: def znodes(self) -> np.ndarray: return self._znodes + # ------------------------------------------------------------------ + # Resize guards: the node arrays are fixed, so the setters refuse + # ------------------------------------------------------------------ + + @property + def origin(self) -> np.ndarray: + """Return the origin of the grid. + + Returns + ------- + np.ndarray + The origin, shape (3,). + """ + return self._origin + + @origin.setter + def origin(self, origin: np.ndarray) -> None: + """Refuse to move the origin. + + Parameters + ---------- + origin : np.ndarray + The requested origin. It is not used. + + Raises + ------ + NotImplementedError + Always. The node arrays are fixed. No state changes. + """ + raise NotImplementedError(_RESIZE_MESSAGE) + + @property + def nsteps(self) -> np.ndarray: + """Return the number of nodes per axis. + + Returns + ------- + np.ndarray + The node counts, shape (3,). + """ + return self._nsteps + + @nsteps.setter + def nsteps(self, nsteps: np.ndarray) -> None: + """Refuse to change the node counts. + + Parameters + ---------- + nsteps : np.ndarray + The requested node counts. It is not used. + + Raises + ------ + NotImplementedError + Always. The node arrays are fixed. No state changes. + """ + raise NotImplementedError(_RESIZE_MESSAGE) + + @property + def step_vector(self) -> np.ndarray: + """Return the step vector of the grid. + + Returns + ------- + np.ndarray + The step vector, shape (3,). + """ + return self._step_vector + + @step_vector.setter + def step_vector(self, step_vector: np.ndarray) -> None: + """Refuse to change the step vector. + + Parameters + ---------- + step_vector : np.ndarray + The requested step vector. It is not used. + + Raises + ------ + NotImplementedError + Always. The node arrays are fixed. No state changes. + """ + raise NotImplementedError(_RESIZE_MESSAGE) + def onGeometryChange(self): if self.interpolator is not None: self.interpolator.reset() diff --git a/packages/loop_common/tests/data/structured_support_reference.npz b/packages/loop_common/tests/data/structured_support_reference.npz index f4330e17..3ecac675 100644 Binary files a/packages/loop_common/tests/data/structured_support_reference.npz and b/packages/loop_common/tests/data/structured_support_reference.npz differ diff --git a/packages/loop_common/tests/test_cell_centre_sampling.py b/packages/loop_common/tests/test_cell_centre_sampling.py new file mode 100644 index 00000000..87dac7e7 --- /dev/null +++ b/packages/loop_common/tests/test_cell_centre_sampling.py @@ -0,0 +1,126 @@ +"""Cell centres of a BoundingBox agree with the cells of its structured grids, also for rotated boxes.""" + +# The task premise for a one-to-one test against the enclosure grid was wrong. +# The evidence for FR-015 and SC-011 is in test_centres_are_cell_midpoints_in_world. +import itertools + +import numpy as np +import pytest +from loop_common.geometry import BoundingBox + +BOXES = { + "unit": dict(origin=[0.0, 0.0, 0.0], maximum=[1.0, 1.0, 1.0], nsteps=[10, 10, 10]), + "offset_anisotropic": dict( + origin=[100.0, -20.0, 5.0], maximum=[112.0, 4.0, 35.0], nsteps=[7, 5, 4] + ), +} +ANGLES = [0.0, 30.0, 90.0] +LOCAL_ORIGINS = {"origin0": (0.0, 0.0, 0.0), "origin_shifted": (7.0, -3.0, 11.0)} + + +def _rotation_about_z(angle): + radians = np.deg2rad(angle) + rotation = np.eye(3) + rotation[:2, :2] = [[np.cos(radians), -np.sin(radians)], [np.sin(radians), np.cos(radians)]] + return rotation + + +@pytest.fixture( + params=list(itertools.product(BOXES, ANGLES, LOCAL_ORIGINS)), + ids=lambda p: "-".join(map(str, p)), +) +def box(request): + name, angle, origin_name = request.param + box = BoundingBox(**BOXES[name]) + box.set_local_transform( + local_origin=np.array(LOCAL_ORIGINS[origin_name]), rotation_matrix=_rotation_about_z(angle) + ) + return box + + +def _tolerance(box): + """Point tolerance, the same as in test_rigid_transform_contract.""" + length = np.linalg.norm(box.maximum - box.origin) + world_magnitude = max(np.linalg.norm(box.origin), np.linalg.norm(box.maximum)) + return 1e-10 * length + 32 * np.finfo(float).eps * max(world_magnitude, length) + + +def _world_centres(box): + return box.reproject(box.cell_centres(order="F")) + + +def test_centres_round_trip_and_index_own_cells(box): + """The centres survive a world round trip and index the cells of the box own grid. + + This is a consistency check on the cell grid of the box itself. The grid + is built from the same linspace as the centres. This test is not evidence + for FR-015 or SC-011. + """ + n_cells = int(np.prod(box.nsteps - 1)) + centres = box.cell_centres(order="F") + # World round trip must give the same local centres. + local = box.project(box.reproject(centres)) + np.testing.assert_allclose(local, centres, atol=_tolerance(box), rtol=0.0) + # The grid nodes are linspace(origin, maximum) in the frame of the box, and + # project() maps them to the local frame. So go back with reproject(), then + # find the cell index on that grid: x changes fastest (order F). + origin = box.regular_grid(order="F", local=False)[0] + step = box.node_step_vector + cell = np.floor((box.reproject(local) - origin) / step).astype(int) + nx, ny, nz = (int(n) - 1 for n in box.nsteps) + assert len(cell) == n_cells + assert np.all(cell >= 0) + assert np.all(cell < np.array([nx, ny, nz])) + index = cell[:, 0] + nx * (cell[:, 1] + ny * cell[:, 2]) + np.testing.assert_array_equal( + np.bincount(index, minlength=n_cells), np.ones(n_cells, dtype=int) + ) + np.testing.assert_array_equal(index, np.arange(n_cells)) + + +def test_local_coordinates_grid_associates_own_centres(box): + """Containment check: each cell centre falls inside the enclosure grid. + + The grid is structured_grid(local_coordinates=True). The check gives no + -1 for any centre. The enclosure grid has axis-aligned cells that differ + from the cells of a rotated box, so a one-to-one count against it does + not hold. + """ + grid = box.structured_grid(local_coordinates=True) + assert np.all(grid.cell_index_at(box.project(_world_centres(box))) != -1) + + +def test_centres_are_cell_midpoints_in_world(box): + """Each world centre is the mean of the 8 corner nodes of its cell. + + This test is evidence for FR-015 and SC-011. Each cell of the exported + grid has exactly one centre, as the mean of its corners. The centres are + all different. + """ + grid = box.structured_grid() + nx, ny, nz = (int(n) for n in box.nsteps) + k, j, i = np.meshgrid(np.arange(nz - 1), np.arange(ny - 1), np.arange(nx - 1), indexing="ij") + i, j, k = i.ravel(), j.ravel(), k.ravel() # x varies fastest, as the flat cell index + corners = np.zeros((len(i), 3)) + for di, dj, dk in itertools.product((0, 1), repeat=3): + corners += grid.nodes[(i + di) + nx * ((j + dj) + ny * (k + dk))] + np.testing.assert_allclose(_world_centres(box), corners / 8.0, atol=_tolerance(box), rtol=0.0) + # No two cells share a centre, so each cell has exactly one centre. + decimals = 6 + assert len(np.unique(np.round(corners / 8.0, decimals), axis=0)) == len(i) + + +def test_centres_match_exported_cells(box): + """The world centres equal the cell centres of the exported grid.""" + grid = box.structured_grid() + np.testing.assert_allclose( + _world_centres(box), grid.cell_centres, atol=_tolerance(box), rtol=0.0 + ) + + +def test_centres_inside_box(box): + """Each world centre is inside the box limits.""" + world = _world_centres(box) + tolerance = _tolerance(box) + assert np.all(world >= box.origin - tolerance) + assert np.all(world <= box.maximum + tolerance) diff --git a/packages/loop_common/tests/test_resize_paths.py b/packages/loop_common/tests/test_resize_paths.py new file mode 100644 index 00000000..6cb3c118 --- /dev/null +++ b/packages/loop_common/tests/test_resize_paths.py @@ -0,0 +1,183 @@ +"""The resize paths of the structured grids keep the extent and the counts. + +Each resize setter must keep the grid extent, or enclose it, and must keep the node counts whole. +Setting the same value again and again must change nothing. The tests give a wrong answer on a +grid that drifts. +""" + +import numpy as np +import pytest +from loop_common.geometry import StructuredGrid2DGeometry, StructuredGrid3DGeometry +from loop_common.supports import RectilinearGrid + +# (origin, nodes per axis, step vector). nsteps counts nodes. +GRIDS_3D = { + "G1": ((0.0, 0.0, 0.0), (11, 6, 4), (1.0, 2.0, 3.0)), + "G2": ((100.0, -20.0, 5.0), (5, 5, 5), (0.5, 0.5, 0.5)), + "G3": ((-1e4, 2e3, -30.0), (21, 3, 8), (10.0, 100.0, 7.5)), +} +GRIDS_2D = { + "H1": ((0.0, 0.0), (5, 5), (1.0, 1.0)), + "H2": ((10.0, -5.0), (8, 4), (0.5, 2.0)), +} +SETTERS = ("origin", "maximum", "nsteps", "step_vector") + + +def _grid_3d(name): + origin, nsteps, step = GRIDS_3D[name] + return StructuredGrid3DGeometry(origin=origin, nsteps=nsteps, step_vector=step) + + +def _state(grid): + """Return a copy of the values that define the grid.""" + return { + "origin": np.array(grid.origin, dtype=float), + "maximum": np.array(grid.maximum, dtype=float), + "step_vector": np.array(grid.step_vector, dtype=float), + "nsteps": np.array(grid.nsteps), + } + + +def _assert_state_equal(grid, expected): + now = _state(grid) + for key, value in expected.items(): + np.testing.assert_array_equal(now[key], value, err_msg=key) + + +def _tolerance(grid_state): + """Length tolerance: 1e-12 * max(1, L) for each axis.""" + length = grid_state["maximum"] - grid_state["origin"] + return 1e-12 * np.maximum(1.0, length) + + +@pytest.mark.parametrize("setter", SETTERS) +@pytest.mark.parametrize("name", list(GRIDS_3D)) +def test_reapplying_same_value_10_times_changes_nothing_3d(name, setter): + grid = _grid_3d(name) + start = _state(grid) + for _ in range(10): + setattr(grid, setter, np.array(start[setter])) + now = _state(grid) + tol = _tolerance(start) + for key in ("origin", "maximum", "step_vector"): + assert np.all(np.abs(now[key] - start[key]) <= tol), (name, setter, key, now[key]) + np.testing.assert_array_equal(now["nsteps"], start["nsteps"]) + + +@pytest.mark.parametrize("name", list(GRIDS_3D)) +def test_nsteps_setter_keeps_extent(name): + grid = _grid_3d(name) + start = _state(grid) + tol = _tolerance(start) + for request in (np.array(start["nsteps"]) * 2, np.array(start["nsteps"])): + grid.nsteps = request + assert np.all(np.abs(grid.maximum - start["maximum"]) <= tol), (name, request, grid.maximum) + assert grid.n_nodes == np.prod(request), (name, request, grid.n_nodes) + + +@pytest.mark.parametrize("name", list(GRIDS_3D)) +def test_step_vector_setter_encloses_extent(name): + grid = _grid_3d(name) + start = _state(grid) + new_step = start["step_vector"] * 1.5 + grid.step_vector = new_step.copy() + tol = _tolerance(start) + assert np.all(grid.maximum >= start["maximum"] - tol), (name, grid.maximum) + assert np.all(grid.maximum < start["maximum"] + new_step), (name, grid.maximum) + + +@pytest.mark.parametrize("name", list(GRIDS_3D)) +def test_step_vector_a_b_a_extent_drift(name): + grid = _grid_3d(name) + start = _state(grid) + step_a = start["step_vector"].copy() + step_b = step_a * 1.5 + grid.step_vector = step_a.copy() + grid.step_vector = step_b.copy() + grid.step_vector = step_a.copy() + drift = grid.maximum - start["maximum"] + print(f"{name}: A-B-A extent drift {drift}, limit {step_a + step_b}") + tol = _tolerance(start) + assert np.all(drift >= -tol), (name, drift) + assert np.all(drift < step_a + step_b), (name, drift) + + +def test_nsteps_below_two_is_rejected_without_change(): + grid = _grid_3d("G1") + start = _state(grid) + with pytest.raises(ValueError): + grid.nsteps = np.array([1, 6, 4]) + _assert_state_equal(grid, start) + + +def test_nsteps_not_whole_is_rejected_without_change(): + grid = _grid_3d("G1") + start = _state(grid) + for bad in ([5.5, 5, 5], [np.nan, 5, 5]): + with pytest.raises(ValueError): + grid.nsteps = np.array(bad) + _assert_state_equal(grid, start) + grid.nsteps = np.array([5.0, 5.0, 5.0]) + assert np.asarray(grid.nsteps).dtype.kind == "i" + np.testing.assert_array_equal(grid.nsteps, [5, 5, 5]) + + +@pytest.mark.parametrize("name", list(GRIDS_3D)) +def test_set_nelements_keeps_extent(name): + grid = _grid_3d(name) + start = _state(grid) + length = start["maximum"] - start["origin"] + h = (np.prod(length) / 200) ** (1 / 3) + result = grid.set_nelements(200) + tol = _tolerance(start) + assert np.all(np.abs(grid.maximum - start["maximum"]) <= tol), (name, grid.maximum) + assert result == grid.n_elements + expected = np.maximum(np.ceil(length / h).astype(int), 2) + np.testing.assert_array_equal(grid.nsteps, expected) + + +def test_set_nelements_small_request_is_clamped(): + grid = _grid_3d("G1") + start = _state(grid) + grid.set_nelements(1) + assert np.all(np.asarray(grid.nsteps) >= 2), grid.nsteps + tol = _tolerance(start) + assert np.all(np.abs(grid.maximum - start["maximum"]) <= tol), grid.maximum + + +@pytest.mark.parametrize("check", ("cells", "n_cell", "n_elements", "dtype", "maximum")) +@pytest.mark.parametrize("name", list(GRIDS_2D)) +def test_2d_counts_are_derived(name, check): + origin, nsteps, step = GRIDS_2D[name] + grid = StructuredGrid2DGeometry(origin=origin, nsteps=nsteps, step_vector=step) + grid.nsteps = [nsteps[0] + 3.0, nsteps[1] + 2.0] + new = np.array([nsteps[0] + 3, nsteps[1] + 2]) + if check == "cells": + np.testing.assert_array_equal(grid.nsteps_cells, new - 1) + elif check == "n_cell": + assert (grid.n_cell_x, grid.n_cell_y) == (new[0] - 1, new[1] - 1) + elif check == "n_elements": + assert grid.n_elements == (new[0] - 1) * (new[1] - 1) + elif check == "dtype": + assert np.asarray(grid.nsteps).dtype.kind == "i" + else: + # The 2D maximum keeps its legacy value: origin + nsteps * step_vector. + np.testing.assert_allclose(grid.maximum, np.array(origin) + new * np.array(step)) + + +@pytest.mark.parametrize("setter", ("nsteps", "step_vector", "origin")) +def test_rectilinear_resize_is_rejected(setter): + grid = RectilinearGrid(np.linspace(0, 5, 6), np.linspace(0, 3, 4), np.linspace(0, 2, 3)) + nodes = np.array(grid.nodes) + origin = np.array(grid.origin) + nsteps = np.array(grid.nsteps) + value = { + "nsteps": np.array([7, 5, 4]), + "step_vector": np.array([2.0, 2.0, 2.0]), + "origin": np.array([1.0, 1.0, 1.0]), + }[setter] + with pytest.raises(NotImplementedError): + setattr(grid, setter, value) + np.testing.assert_array_equal(grid.nodes, nodes) + np.testing.assert_array_equal(grid.origin, origin) + np.testing.assert_array_equal(grid.nsteps, nsteps) diff --git a/packages/loop_interpolation/src/loop_interpolation/__init__.py b/packages/loop_interpolation/src/loop_interpolation/__init__.py index fac883f5..608bef59 100644 --- a/packages/loop_interpolation/src/loop_interpolation/__init__.py +++ b/packages/loop_interpolation/src/loop_interpolation/__init__.py @@ -13,6 +13,7 @@ "DirectionalRegularisation", "DiscreteFoldInterpolator", "DiscreteInterpolator", + "ExclusionRecord", "FDFoldInterpolator", "FiniteDifferenceInterpolator", "FoldEvent", @@ -31,12 +32,15 @@ "PiecewiseLinearInterpolator", "RegionCoverageDiagnostics", "RegularisationConfig", + "ResidualSummary", "StaleFrameError", "StructuredGrid", "StructuredGrid2D", "SurfeRBFInterpolator", "TetMesh", "UnStructuredTetMesh", + "build_exclusion_records", + "direction_cosines", "get_fold_rotation_profile", ] from loop_common.logging import get_logger as getLogger @@ -60,7 +64,11 @@ from ._diagnostics import ( ConstraintDiagnosticsReport, ConstraintFamilyDiagnostics, + ExclusionRecord, RegionCoverageDiagnostics, + ResidualSummary, + build_exclusion_records, + direction_cosines, ) from ._discrete_fold_interpolator import ( DiscreteFoldInterpolator, diff --git a/packages/loop_interpolation/src/loop_interpolation/_diagnostics.py b/packages/loop_interpolation/src/loop_interpolation/_diagnostics.py index f1f591ce..8481639d 100644 --- a/packages/loop_interpolation/src/loop_interpolation/_diagnostics.py +++ b/packages/loop_interpolation/src/loop_interpolation/_diagnostics.py @@ -1,6 +1,17 @@ from __future__ import annotations -from dataclasses import dataclass, field +from dataclasses import asdict, dataclass, field + +import numpy as np + +EXCLUSION_EXPLANATIONS: dict[str, str] = { + "invalid_reference_value": "The reference value of the point is not a finite number.", + "invalid_reference_direction": "The reference direction of the point is zero or not finite.", + "evaluation_error": "The model could not be evaluated at the point.", + "outside_model": "The point is outside the model bounds. The count is in points.", + "unattributed": "The rows were dropped for an unknown reason. The count is in matrix rows.", +} +UNKNOWN_EXPLANATION = "No explanation is registered for this reason." @dataclass(frozen=True) @@ -24,12 +35,185 @@ class RegionCoverageDiagnostics: active_fraction: float +@dataclass(frozen=True) +class ResidualSummary: + """Statistics of the residuals of one constraint family. + + Attributes + ---------- + name : str + Name of the residual set. + kind : str + Either ``"physical"`` or ``"assembled"``. + unit : str + One of ``"field"``, ``"degree"`` or ``"assembled"``. + count : int + Number of finite residuals. + nan_count : int + Number of residuals that are not finite. + rms : float or None + Root mean square of the finite residuals. None when ``count`` is 0. + max_abs : float or None + Largest absolute finite residual. None when ``count`` is 0. + mean_abs : float or None + Mean absolute finite residual. None when ``count`` is 0. + """ + + name: str + kind: str + unit: str + count: int + nan_count: int + rms: float | None + max_abs: float | None + mean_abs: float | None + + +@dataclass(frozen=True) +class ExclusionRecord: + """One reason why rows or points were excluded. + + Attributes + ---------- + family : str + Name of the constraint family. + reason : str + Short key of the reason. + count : int + Number of excluded items. + explanation : str + One sentence that explains the reason. + """ + + family: str + reason: str + count: int + explanation: str + + +def summarise_residuals(name: str, kind: str, unit: str, residuals: np.ndarray) -> ResidualSummary: + """Summarise a set of residuals. + + Parameters + ---------- + name : str + Name of the residual set. + kind : str + Either ``"physical"`` or ``"assembled"``. + unit : str + One of ``"field"``, ``"degree"`` or ``"assembled"``. + residuals : numpy.ndarray + The residuals. The array is flattened and is not changed. + + Returns + ------- + ResidualSummary + Statistics of the finite residuals. Values that are not finite + are only counted. + """ + flat = np.asarray(residuals, dtype=float).ravel() + finite = flat[np.isfinite(flat)] + nan_count = int(flat.size - finite.size) + if finite.size == 0: + return ResidualSummary(name, kind, unit, 0, nan_count, None, None, None) + absolute = np.abs(finite) + return ResidualSummary( + name=name, + kind=kind, + unit=unit, + count=int(finite.size), + nan_count=nan_count, + rms=float(np.sqrt(np.mean(finite**2))), + max_abs=float(np.max(absolute)), + mean_abs=float(np.mean(absolute)), + ) + + +def direction_cosines(model_vectors: np.ndarray, data_vectors: np.ndarray) -> np.ndarray: + """Compute the cosine of the angle between paired vectors. + + Parameters + ---------- + model_vectors : numpy.ndarray + Array of shape (n, d). + data_vectors : numpy.ndarray + Array of shape (n, d). + + Returns + ------- + numpy.ndarray + Array of shape (n,). A row with a zero or not finite norm gives NaN. + + Raises + ------ + ValueError + If the two arrays do not have the same shape. + """ + model = np.asarray(model_vectors, dtype=float) + data = np.asarray(data_vectors, dtype=float) + if model.shape != data.shape: + raise ValueError(f"Shape mismatch: model vectors {model.shape}, data vectors {data.shape}.") + model_norm = np.linalg.norm(model, axis=-1) + data_norm = np.linalg.norm(data, axis=-1) + valid = np.isfinite(model_norm) & np.isfinite(data_norm) & (model_norm > 0) & (data_norm > 0) + result = np.full(model_norm.shape, np.nan) + dot = np.sum(model * data, axis=-1) + result[valid] = dot[valid] / (model_norm[valid] * data_norm[valid]) + return result + + +def build_exclusion_records( + excluded_rows: dict[str, dict[str, int]], + families: dict[str, ConstraintFamilyDiagnostics], + rows_per_point: dict[str, int] | None = None, +) -> tuple[ExclusionRecord, ...]: + """Build the records of excluded rows and points. + + Parameters + ---------- + excluded_rows : dict + Maps a family name to a dict of reason and count. + families : dict + Maps a family name to its diagnostics. + rows_per_point : dict, optional + Number of matrix rows that one point makes, for each family. + The default is 1 for each family. + + Returns + ------- + tuple of ExclusionRecord + Records sorted by family and reason. A count of 0 makes no record. + """ + rows_per_point = rows_per_point or {} + records = [] + + def add(family: str, reason: str, count: int) -> None: + """Add one record when the count is more than 0.""" + if count > 0: + explanation = EXCLUSION_EXPLANATIONS.get(reason, UNKNOWN_EXPLANATION) + records.append(ExclusionRecord(family, reason, int(count), explanation)) + + for family, reasons in excluded_rows.items(): + for reason, count in reasons.items(): + add(family, reason, count) + for name, diagnostics in families.items(): + outside = diagnostics.outside_model_point_count + add(name, "outside_model", outside) + if diagnostics.dropped_rows is not None: + k = rows_per_point.get(name, 1) + add(name, "unattributed", diagnostics.dropped_rows - outside * k) + return tuple(sorted(records, key=lambda r: (r.family, r.reason))) + + @dataclass(frozen=True) class ConstraintDiagnosticsReport: interpolator_type: str families: dict[str, ConstraintFamilyDiagnostics] = field(default_factory=dict) region_coverage: RegionCoverageDiagnostics | None = None outside_model_points: dict[str, int] = field(default_factory=dict) + physical_residuals: dict[str, ResidualSummary] = field(default_factory=dict) + assembled_residuals: dict[str, ResidualSummary] = field(default_factory=dict) + exclusions: tuple[ExclusionRecord, ...] = () @property def total_rows(self) -> int: @@ -40,7 +224,7 @@ def active_families(self) -> dict[str, ConstraintFamilyDiagnostics]: return {k: v for k, v in self.families.items() if v.active} def to_dict(self) -> dict: - return { + data = { "interpolator_type": self.interpolator_type, "total_rows": self.total_rows, "families": { @@ -56,16 +240,29 @@ def to_dict(self) -> dict: } for name, family in self.families.items() }, - "region_coverage": None - if self.region_coverage is None - else { - "total_support_nodes": self.region_coverage.total_support_nodes, - "active_region_nodes": self.region_coverage.active_region_nodes, - "inactive_region_nodes": self.region_coverage.inactive_region_nodes, - "active_fraction": self.region_coverage.active_fraction, - }, + "region_coverage": ( + None + if self.region_coverage is None + else { + "total_support_nodes": self.region_coverage.total_support_nodes, + "active_region_nodes": self.region_coverage.active_region_nodes, + "inactive_region_nodes": self.region_coverage.inactive_region_nodes, + "active_fraction": self.region_coverage.active_fraction, + } + ), "outside_model_points": dict(self.outside_model_points), } + if self.physical_residuals: + data["physical_residuals"] = { + key: asdict(value) for key, value in self.physical_residuals.items() + } + if self.assembled_residuals: + data["assembled_residuals"] = { + key: asdict(value) for key, value in self.assembled_residuals.items() + } + if self.exclusions: + data["exclusions"] = [asdict(record) for record in self.exclusions] + return data def summary(self) -> str: lines = [ @@ -99,4 +296,20 @@ def summary(self) -> str: for key, value in sorted(self.outside_model_points.items()): lines.append(f" - {key}: {value}") + for title, block in ( + ("Physical residuals:", self.physical_residuals), + ("Assembled residuals:", self.assembled_residuals), + ): + if block: + lines.append(title) + for key, item in sorted(block.items()): + lines.append( + f" - {key}: count={item.count}, nan={item.nan_count}, " + f"rms={item.rms}, max_abs={item.max_abs} ({item.unit})" + ) + if self.exclusions: + lines.append("Exclusions:") + for record in self.exclusions: + lines.append(f" - {record.family}: {record.reason}={record.count}") + return "\n".join(lines) diff --git a/packages/loop_interpolation/tests/test_residual_exclusion_report.py b/packages/loop_interpolation/tests/test_residual_exclusion_report.py new file mode 100644 index 00000000..0985a117 --- /dev/null +++ b/packages/loop_interpolation/tests/test_residual_exclusion_report.py @@ -0,0 +1,233 @@ +import dataclasses + +import numpy as np +import pytest +from loop_interpolation._diagnostics import ( + UNKNOWN_EXPLANATION, + ConstraintDiagnosticsReport, + ConstraintFamilyDiagnostics, + ExclusionRecord, + RegionCoverageDiagnostics, + ResidualSummary, + build_exclusion_records, + direction_cosines, + summarise_residuals, +) + + +def _family( + name, + dropped_rows=0, + source_point_count=0, + outside_model_point_count=0, + row_count=10, +): + return ConstraintFamilyDiagnostics( + name=name, + active=True, + row_count=row_count, + dropped_rows=dropped_rows, + effective_weight_mean=1.0, + effective_weight_min=1.0, + effective_weight_max=1.0, + source_point_count=source_point_count, + outside_model_point_count=outside_model_point_count, + ) + + +def _unattributed(records): + return [r for r in records if r.reason == "unattributed"] + + +def test_direction_cosines_pure_function(): + model = np.array( + [ + [1.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [0.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + ] + ) + data = np.array( + [ + [2.0, 0.0, 0.0], + [-3.0, 0.0, 0.0], + [0.0, 4.0, 0.0], + [1.0, 0.0, 0.0], + [np.nan, 0.0, 0.0], + ] + ) + model_before = model.copy() + data_before = data.copy() + result = direction_cosines(model, data) + assert result.shape == (5,) + np.testing.assert_allclose(result[:3], [1.0, -1.0, 0.0], atol=1e-12) + assert np.isnan(result[3]) + assert np.isnan(result[4]) + # Inputs are not changed. + np.testing.assert_array_equal(model, model_before) + np.testing.assert_array_equal(data, data_before) + # A shape mismatch is an error. + with pytest.raises(ValueError): + direction_cosines(model[:3], data[:2]) + with pytest.raises(ValueError): + direction_cosines(model, data[:, :2]) + # The result does not depend on the length of the vectors. + a = np.array([[1.0, 2.0, 3.0], [0.0, 1.0, 1.0]]) + b = np.array([[3.0, 1.0, -2.0], [1.0, 1.0, 0.0]]) + np.testing.assert_allclose( + direction_cosines(a, b), direction_cosines(5.0 * a, 0.25 * b), atol=1e-12 + ) + + +def test_report_without_verification_is_unchanged(): + family = _family("value", source_point_count=4) + region = RegionCoverageDiagnostics( + total_support_nodes=10, + active_region_nodes=6, + inactive_region_nodes=4, + active_fraction=0.6, + ) + report = ConstraintDiagnosticsReport( + interpolator_type="FiniteDifferenceInterpolator", + families={"value": family}, + region_coverage=region, + outside_model_points={"value": 1}, + ) + data = report.to_dict() + assert set(data.keys()) == { + "interpolator_type", + "total_rows", + "families", + "region_coverage", + "outside_model_points", + } + assert data["total_rows"] == 10 # row_count of the one family + assert set(data["families"]["value"].keys()) == { + "active", + "row_count", + "dropped_rows", + "effective_weight_mean", + "effective_weight_min", + "effective_weight_max", + "source_point_count", + "outside_model_point_count", + } + assert set(data["region_coverage"].keys()) == { + "total_support_nodes", + "active_region_nodes", + "inactive_region_nodes", + "active_fraction", + } + # The old fields keep their order. New fields can only come after them. + assert [f.name for f in dataclasses.fields(ConstraintDiagnosticsReport)][:4] == [ + "interpolator_type", + "families", + "region_coverage", + "outside_model_points", + ] + assert [f.name for f in dataclasses.fields(ConstraintFamilyDiagnostics)][:9] == [ + "name", + "active", + "row_count", + "dropped_rows", + "effective_weight_mean", + "effective_weight_min", + "effective_weight_max", + "source_point_count", + "outside_model_point_count", + ] + + +def test_builder_records_keep_their_own_key(): + # The builder key is not a family that counts it in dropped_rows. + # So dropped_rows is 0 and the rule gives no unattributed record. + families = {"restricted_orthogonal": _family("restricted_orthogonal", dropped_rows=0)} + before = dataclasses.replace(families["restricted_orthogonal"]) + records = build_exclusion_records( + {"restricted_orthogonal": {"invalid_reference_direction": 3}}, families + ) + assert len(records) == 1 + record = records[0] + assert isinstance(record, ExclusionRecord) + assert record.family == "restricted_orthogonal" + assert record.reason == "invalid_reference_direction" + assert record.count == 3 + assert _unattributed(records) == [] + assert families["restricted_orthogonal"] == before + + +def test_unknown_reason_keeps_its_count(): + # dropped_rows is 0, so the rule gives no unattributed record. + families = {"value": _family("value", dropped_rows=0)} + records = build_exclusion_records({"value": {"x": 2}}, families) + assert len(records) == 1 + assert records[0].reason == "x" + assert records[0].count == 2 + assert records[0].explanation == UNKNOWN_EXPLANATION + + +def test_unattributed_rows_are_reported(): + families = { + "value": _family( + "value", dropped_rows=5, source_point_count=10, outside_model_point_count=2 + ) + } + records = build_exclusion_records({}, families) + unattributed = _unattributed(records) + assert len(unattributed) == 1 + assert unattributed[0].family == "value" + assert unattributed[0].count == 3 + + +def test_gradient_outside_points_give_no_false_unattributed(): + def make(dropped): + return { + "gradient": _family( + "gradient", + dropped_rows=dropped, + source_point_count=10, + outside_model_point_count=2, + ) + } + + records = build_exclusion_records({}, make(4), rows_per_point={"gradient": 2}) + # 4 dropped rows = 2 outside points * 2 rows. No unattributed record. + # One outside_model record gives the count in points. + assert _unattributed(records) == [] + assert len(records) == 1 + assert records[0].family == "gradient" + assert records[0].reason == "outside_model" + assert records[0].count == 2 + records = build_exclusion_records({}, make(6), rows_per_point={"gradient": 2}) + unattributed = _unattributed(records) + assert len(unattributed) == 1 + assert unattributed[0].count == 2 + + +def test_dropped_rows_none_gives_no_unattributed(): + families = {"interface": _family("interface", dropped_rows=None)} + assert build_exclusion_records({}, families) == () + + +def test_surfe_style_dropped_equals_outside(): + families = { + "value": _family( + "value", dropped_rows=2, source_point_count=10, outside_model_point_count=2 + ) + } + records = build_exclusion_records({}, families, rows_per_point=None) + # Dropped rows equal outside points. No unattributed record. + # One outside_model record gives the count. + assert _unattributed(records) == [] + assert len(records) == 1 + assert records[0].family == "value" + assert records[0].reason == "outside_model" + assert records[0].count == 2 + + +def test_residual_summary_is_importable(): + # The name must exist. Its behaviour is tested with the implementation task. + assert ResidualSummary is not None + assert callable(summarise_residuals) diff --git a/specs/005-coordinate-frame-correctness/plan.md b/specs/005-coordinate-frame-correctness/plan.md index 41b59561..44776f1f 100644 --- a/specs/005-coordinate-frame-correctness/plan.md +++ b/specs/005-coordinate-frame-correctness/plan.md @@ -7,18 +7,20 @@ The constitution (.specify/memory/constitution.md) has priority. --> ## Open points -- **OP1**: This plan gives the exact design only for stages S1 to S4 (S1: +- **OP1**: This plan gives the exact design only for stages S1 to S5 (S1: shared transform, reference-constraint boundary and fault frame; S2: fold, P2, directional-regularisation and constant-norm cells; S3: saved frames, copies and transform changes; S4: display, data export, fault-centre geometry and - intrusion classification). The spec-planner limit is 10 code files, and the + intrusion classification; S5: resize paths, cell centres, non-mutating + diagnostics, residual and exclusion report). The spec-planner limit is 10 code files, and the maintainer selected a staged plan on 2026-11-04. S2 has its exact design since 2026-11-04 ("Stage S2 design"). S3 has its exact design since 2026-11-05 ("Stage S3 design"), and S4 since 2026-11-05 ("Stage S4 design"). - Stages S5 and S6 have a stage map, a file list from a grep survey and a - requirement allocation, but no exact signatures. Owner: the spec-planner - skill, one run for each stage. Blocks: the tasks of S5 and S6. The spec-tasks - skill makes the S4 tasks next (S1 is merged; S2 is in pull request #396; S3 + S5 has its exact design since 2026-11-05 ("Stage S5 design"). Stage S6 has + a stage map, a file list from a grep survey and a requirement allocation, but + no exact signatures. Owner: the spec-planner + skill, one run for each stage. Blocks: the tasks of S6. The spec-tasks + skill makes the S5 tasks next (S1 is merged; S2 is in pull request #396; S3 is built in pull request #397). - **OP2**: Run time of the mandatory R2 and R5 matrices. R2 needs 4 anchors x 3 rotations x 2 centre modes x 2 engines = 48 solves (FDI with 20000 @@ -97,6 +99,40 @@ The constitution (.specify/memory/constitution.md) has priority. --> defect. Owner: the maintainer. Decision needed: file an issue, as for the C-005 and C-006 defects. Default until an answer: S4 keeps the expression and tests only the branch without `value_map`. Blocks: nothing. +- **OP12**: The spec (FR-015) says that cell-centre sampling can place a centre + beyond its cell. By reading, `BoundingBox.cell_centres` is correct for a + rotated box. The spacing equals `node_step_vector`, and the four callers + convert with `reproject`. The planner did not run the code. Owner: the S5 + cell-centre test task. Decision needed: if all 60 cells of the new tests pass, + record "not reproduced" and keep the tests as guards. Default until an answer: + no code change in `_bounding_box.py`. If a cell fails, use the closed form of + "S5 decisions" item 7. Blocks: nothing. +- **OP13** (closed 2026-11-05, maintainer answer Q3, see "Answers"): The misfit + methods must not change state and must not rebuild a stale model. Rule: on a + stale builder, fault builder or interpolator, the methods write one + `UserWarning` and return NaN ("S5 decisions" item 11). +- **OP14** (closed 2026-11-05, maintainer answer Q1, see "Answers"): The 2D + `maximum` keeps its legacy value `origin + nsteps * step`. It counts nodes as + cells. S5 keeps it for FR-028. The four existing assertions pass with no change. + S5 only derives `nsteps_cells`, `n_cell_x` and `n_cell_y` from `nsteps`. +- **OP15** (closed 2026-11-05, maintainer answer Q2, see "Answers"): The old misfit + methods keep their names and public shape. New names give the 1-D result. The + gradient value changes meaning. `COMPAT.md` and `ROADMAP.md` record it. +- **OP16**: The old `evaluate_gradient_misfit` for two families (gradient and + normal) of different length. Before S5 it raised `ValueError` on current NumPy. + Equal lengths gave shape `(2, N)`. Owner: the maintainer. Decision needed: keep + `(2, N)` for equal lengths and return one row `(1, N1 + N2)` for different + lengths, or always return one row. Default until an answer: the first choice, + because it keeps the shape that worked. `COMPAT.md` records it. Blocks: + nothing. +- **OP17**: The setter fix and the extent from `set_nelements` can change numeric + tests, including `packages/loop_interpolation/tests/golden` (the parity oracle + for the Rust path, with its own `--update-golden` option). Owner: the S5 survey + task. Decision needed: whether to regenerate such a file. Default until an + answer: record each change in "S5 golden review (C-004)"; regenerate a file only + when the cause is the intended extent change, and name the file; leave the + gates red and report if the cause is not clear. Blocks: the merge of S5 when a + parity file changes. ## Answers @@ -185,6 +221,66 @@ The constitution (.specify/memory/constitution.md) has priority. --> fail before T032 and pass after it. T030 and T031 also need T032. The PLI control meets SC-002 (value error 1.8e-12, gradient error 1.3e-10). +- 2026-11-05 (S5 plan): OP1 is closed for S5. After an independent review of + the S5 design, the maintainer answered five questions on 2026-11-05. The S5 + design was revised to match. +- 2026-11-05 (Q1, 2D `maximum`, maintainer): do not change the value of the 2D + `maximum`. Derive `nsteps_cells`, `n_cell_x` and `n_cell_y` from `nsteps`, and + keep `maximum = origin + nsteps * step` (legacy; it counts nodes as cells). This + is kept for FR-028. The four existing assertions in + `packages/loop_common/tests/test_2d_discrete_support.py` (lines 14 and 22) and + `tests/unit/interpolator/test_2d_discrete_support.py` (lines 16 and 24) stay + unchanged and must pass. OP14 is closed. +- 2026-11-05 (Q2, misfit result, maintainer): keep the old public names and the + old public shape of `evaluate_value_misfit` and `evaluate_gradient_misfit`. The + new non-mutating, world-correct 1-D result has new names + (`evaluate_gradient_cosines`, `evaluate_value_residuals`). The old methods call + the same corrected code. The gradient value changes meaning (a cosine). This is + a behaviour change, recorded in `COMPAT.md` (migration notice) and in a + `ROADMAP.md` decision entry. OP15 is closed. +- 2026-11-05 (Q3, stale model, maintainer): the misfit methods change no state + and do not rebuild a stale model. S5 chose one rule: warn once with + `UserWarning` and return NaN. OP13 is closed. +- 2026-11-05 (Q4, `set_nelements`, maintainer): accept the exact extent after + `set_nelements`. S5 adds a numeric-test survey, a new section "S5 golden review + (C-004)", and a regeneration of the golden data files after review (C-004, + SC-015). The earlier sentence that the 3D golden test passes with no change is + withdrawn. `set_nelements` clamps `nsteps` to at least 2; the setter rejects + non-whole values; the 2D `nsteps` assignment coerces to an int array. +- 2026-11-05 (Q5, value misfit, maintainer): the public value misfit keeps its + normalisation by the field range. The physical residual uses an unnormalised + private helper, so its unit "field" is true. +- 2026-11-05 (S5 build start, orchestrator, five gaps of the task texts, settled + by reading the code): (1) `GeologicalFeature` and its base class have no + `bounding_box` attribute. The misfit code uses `self.interpolator.bounding_box` + (an attribute of `GeologicalInterpolator`, `None` when no frame is set), as the + old code does. `self.model.bounding_box` is not used: the model can be `None`. + (2) A stale model for T081: for a stale interpolator, assign + `interpolator.up_to_date = False` after the solve (a plain attribute). For a + stale builder, assign `builder._up_to_date = False` after the solve. Do not use + `builder.set_not_up_to_date`: the `GeologicalFeatureBuilder` override also sets + the interpolator stale. The test asserts the flag value before the call. + (3) `n_elements` exists on `StructuredGrid2DGeometry` (`nsteps_cells[0] * + nsteps_cells[1]`, `_structured_grid_2d.py:45`) and on the 3D geometry (line + 166). (4) No existing test checks the keys of `ConstraintDiagnosticsReport.to_dict`. + T071 writes the pre-S5 key set as a literal from the code: report keys + `interpolator_type`, `total_rows`, `families`, `region_coverage`, + `outside_model_points`; family keys `active`, `row_count`, `dropped_rows`, + `effective_weight_mean`, `effective_weight_min`, `effective_weight_max`, + `source_point_count`, `outside_model_point_count`; region keys + `total_support_nodes`, `active_region_nodes`, `inactive_region_nodes`, + `active_fraction`; report field order `interpolator_type`, `families`, + `region_coverage`, `outside_model_points`. (5) The gradient family adds 2 + rows for each inside point for FDI (strike and dip rows, `name="gradient"`). + PLI is different: `P1Interpolator.add_gradient_constraints` is `pass` + (`_p1interpolator.py:70`), so PLI assembles no `gradient` rows, and the builder + reports `dropped_rows = 2 * source_points` for the `gradient` family of PLI. + So `verification_report` gives an `unattributed` record for PLI gradient data. + That record is true (PLI ignores gradient data) and it is not a defect of S5. + T081 checks the no-false-unattributed case with FDI only, and for PLI it records + the measured `unattributed` count (expected `2 * source_points`) and the + `outside_model` count. The worker reports both numbers. + ## Summary Observations, fault geometry, masks and public queries stay in world @@ -207,8 +303,9 @@ S1 repairs the shared boundary and the fault frame: S2 repairs the fold paths. S3 makes saved interpolators, copies and transform changes keep the frame ("Stage S3 design"). S4 removes two extra conversions: the display and export of observation data, and the intrusion classification. It also -deprecates `GeologicalModel.rescale` ("Stage S4 design"). S5 and S6 repair grid and diagnostics, and the auxiliary -utilities and validation matrix (see "Stages and rollback"). +deprecates `GeologicalModel.rescale` ("Stage S4 design"). S5 repairs the resize paths and the misfit diagnostics and adds the residual and +exclusion report ("Stage S5 design"). S6 covers the auxiliary utilities and the +validation matrix (see "Stages and rollback"). ## Technical context @@ -274,7 +371,7 @@ utilities and validation matrix (see "Stages and rollback"). ## Design -S1 to S4 have the exact design below. S5 and S6 have scope only (OP1). +S1 to S5 have the exact design below. S6 has scope only (OP1). ### Types and data @@ -1686,6 +1783,862 @@ two codes give the same `get_data()` output; (d) on the S3 code `model.rescale(p and `model.scale(p)` give no warning, on the S4 code each gives one, and the values are equal. +### Stage S5 design + +S5 has four aims. First, S5 makes the remaining resize paths keep the node counts +and cell counts consistent (FR-014). Second, S5 checks that cell-centre sampling +agrees with the grid cells (FR-015). Third, S5 makes the feature misfit +diagnostics leave their inputs unchanged and compare quantities in one basis +(FR-023). Fourth, S5 adds a residual and exclusion report that keeps physical +residuals and assembled residuals apart (FR-024). + +Base: branch `feat/005-s4-world-locations` (S4, pull request #401, not merged). If +S4 changes before it merges, run this planner again on the merged result. + +The plan was revised on 2026-11-05 after an independent review and five +maintainer answers (see "Answers"). The text below is the revised design. + +#### S5 facts from the code (read on 2026-11-05) + +- Resize, 3D. S1 repaired the `origin` and `maximum` setters of + `StructuredGrid3DGeometry` with `_cells_for_length`. Two setters of the same + class were not repaired (`_structured_grid_3d.py`, `nsteps` setter at line 71, + `step_vector` setter at line 131). The `nsteps` setter counts nodes. It does + `step_vector /= nsteps_new / nsteps_old`. This keeps the product + `nsteps * step` constant. It does not keep the extent `(nsteps - 1) * step`. + Example: the node count 11 to 21 with step 1 gives step `0.524` and a new + extent `10.48` instead of `10`. +- The `step_vector` setter sets the node count to `ceil(nsteps / factor)`. It + treats the node count as a cell count. So the extent also drifts. It can lose + a cell. +- `set_nelements` (`_structured_grid_3d.py:50-57`) ends with `self.nsteps = + nsteps`. It computes `nsteps = ceil((maximum - origin) / step)` with a cubic + step from the requested element count. So it uses the `nsteps` setter, and the + setter fix changes what `set_nelements` leaves behind. Callers: + `GeologicalFeatureBuilder.build` calls `interpolator.set_nelements` (line 575 + of `_geological_feature_builder.py`) and the support method is reached through + `DiscreteInterpolator.set_nelements` (line 104). `set_interpolation_geometry` + runs `support.step_vector = support.step_vector * 0.9` in a loop while + `interpolator.dof < 100` (line 525 of the same file). So both setters run in + ordinary model builds, and numeric tests that depend on a built support can + change. +- The golden test `packages/loop_common/tests/test_structured_support_reference.py` + runs `set_nelements(200)`, then `nsteps + 1`, then `step_vector * 2`, then + `origin - 1`, then `maximum + 3`, in sequence on one grid (lines 171 to 180). + Each case records the state that the previous setters left. So the cases + `after_set_nelements`, `after_nsteps`, `after_step_vector`, `after_origin` and + `after_maximum` all change when the two setters change. The earlier text of + this plan said that `after_origin` and `after_maximum` stay: that was wrong. +- `StructuredGrid2DGeometry` (`_structured_grid_2d.py`) has plain attributes + `nsteps`, `maximum`, `nsteps_cells`, `n_cell_x` and `n_cell_y`, set once in + `__init__` (lines 23 to 30). `maximum = origin + nsteps * step_vector` counts + nodes as cells. `nodes` uses `origin + nsteps_cells * step_vector`. The two + disagree by one step. The three count attributes and `maximum` go stale if + `nsteps` is assigned later. The 3D class computes `maximum` as + `origin + nsteps_cells * step_vector`. The 2D class `set_nelements` raises + `NotImplementedError` (`supports/_2d_structured_grid.py:50`). Four existing + assertions read the 2D `maximum`: `packages/loop_common/tests/test_2d_discrete_support.py` + lines 14 and 22, and `tests/unit/interpolator/test_2d_discrete_support.py` + lines 16 and 24. They expect `origin + nsteps * step` (the legacy value). +- `RectilinearGrid` (`supports/_3d_rectilinear_grid.py`) holds node arrays + `_xnodes`, `_ynodes`, `_znodes`. It inherits the `origin`, `nsteps` and + `step_vector` setters of the 3D geometry. A call changes the stored counts and + the stored origin, but not the node arrays. So the geometry and the nodes + disagree. `maximum` is read from the node arrays and has no setter. The + constructors of the geometry classes write the private attributes (`_origin`, + `_nsteps`, `_step_vector`), not the setters. +- Cell centres. `BoundingBox.cell_centres` (`_bounding_box.py:647`) builds + `nsteps - 1` centres with `linspace(origin + h, maximum - h, nsteps - 1)` and + `h = node_step_vector / 2`. Now `maximum - origin = (nsteps - 1) * step`. So + the spacing is exactly `node_step_vector`. Every centre is the middle of a cell + of `regular_grid(local=False)`. The method returns local points. All four + callers (`_base_geological_feature.py:391`, `420`, `442` and + `_model_exporter.py:77`) call `reproject` on the result. The existing tests + (`test_bounding_box_structured_grid.py`, `test_scalar_field_grid.py`) use + unrotated boxes only. By reading, the default path has no defect for a rotated + box (OP12). `structured_grid(local_coordinates=True)` builds a different grid. + Its cells are the enclosure of the projected corners. That grid works in the + local frame, so `grid.cell_index_at` takes local points. Its `cell_centres` are + the centres of those cells. They are not `box.cell_centres()`. The supports' + `cell_centres(global_index)` (3D, 2D, rectilinear) use their own local frame and + do not apply `rotation_xy`, as their docstrings state. +- Misfit methods (`_geological_feature.py`: `evaluate_gradient_misfit` at line 240 + and `evaluate_value_misfit` at line 272 on this branch). `get_gradient_constraints()` + and `get_norm_constraints()` return the stored arrays by reference + (`_geological_interpolator.py`: the definitions are at lines 633 and 652, and + the returns are at lines 640 and 660). So `grad /= norm` (line 255) changes + `interpolator.data["gradient"]` in place. The norm is taken over all columns + (position, vector and weight). So the position and the weight columns change + too. `grad[:, :3:6]` (line 260) is the slice with start 0, stop 3 and step 6. + That is column `0` only. The vector columns are `3:6`. The stored vectors are + in the local basis. The interpolator converts them once when it stores them. + `evaluate_gradient` returns world vectors. So a rotated frame compares two + bases. +- Evaluation path. `GeologicalFeature.evaluate_value` (line 117) and + `evaluate_gradient` call `self.builder.up_to_date()` (lines 141 and 191). + `BaseBuilder.up_to_date` (`_base_builder.py:81`) first calls itself on every + fault builder. Then, if the builder or the interpolator is stale, it calls + `update()`, which builds again and solves. So a misfit method that calls + `self.evaluate_value` on a stale model rebuilds it. This also holds for + `self.max()` and `self.min()` (`_base_geological_feature.py:280` and `294`), + because they call `evaluate_value`. The earlier text of this plan marked this + path as not read. It is now read. No separate cache-validity flag exists in + `GeologicalFeature`, the builder or the interpolator (grep for `cache` on + 2026-11-05). The state is `interpolator.up_to_date`, `builder._up_to_date` + (a protected attribute; no public property exists), `interpolator.c` and + `interpolator.latest_diagnostics_report`. +- Result shape. With one family, `np.array([list])` gives the shape `(1, N)`. + With two families of equal length it gives `(2, N)`. With two families of + different length, `np.array([list, ndarray])` raises `ValueError` on current + NumPy. The value misfit divides by `self.max() - self.min()` (line 288): zero + range gives a division by zero. It does not change its inputs. No code in + `LoopStructural` or `packages` calls the two methods (`git grep -n -w` on + 2026-11-05). `ROADMAP.md` line 832 names them in a history note only. +- Diagnostics report. `_diagnostics.py` holds three frozen dataclasses and + `ConstraintDiagnosticsReport` with `to_dict` and `summary`. It has + `from __future__ import annotations`. The two builders + (`GeologicalInterpolator._build_constraint_diagnostics_report`, line 731, and the + `DiscreteInterpolator` override, line 221) only read state. + `get_constraint_diagnostics_report(refresh=False)` stores the result in + `latest_diagnostics_report`. +- `DiscreteInterpolator.calculate_residual_for_constraints()` (line 497) gives + `A x - b` for each assembled constraint, as a dict by constraint name. The names + are the constraint names: FDI uses `value` and `gradient`, PLI uses `value` and + `norm`. It does not test `self.c`. The constructor sets `self.c` to zeros + (lines 64 to 68), so `c` is not `None` after construction. Line 1355 can + assign another value. +- Dropped rows. The `DiscreteInterpolator` builder sets `dropped_rows` as + `source_points - rows` for `value`, `normal`, `tangent` and `inequality_value` + (line 265), and as `2 * source_points - rows` for `gradient` (line 267). Other + families (interface, inequality pairs, regularisation) get `None`. So + `dropped_rows` is in matrix rows. `outside_model_point_count` is in points. The + base builder (used by `SurfeRBFInterpolator`, which does not override it) sets + `dropped_rows = outside` and `row_count = point count`, so for Surfe both are in + points. +- Builder exclusions. `GeologicalFeatureBuilder.excluded_rows` + (`dict[str, dict[str, int]]`, S1, line 88) uses constraint names as keys: + `_orthogonal` (line 343) and `_equality` (lines 387 + and 393). The reasons are `invalid_reference_value`, + `invalid_reference_direction` and `evaluation_error`. The S1 restricted-support + case uses the feature name `restricted`, so the keys are `restricted_equality` + and `restricted_orthogonal`. These keys are not family names. No report reads + them. +- Python version. `pyproject.toml` and `packages/loop_common/pyproject.toml` set + `requires-python >=3.9`. `_structured_grid_3d.py` has no `from __future__ + import annotations`. The syntax `np.ndarray | None` fails there on Python 3.9. +- Documents. `git grep` finds no mention of `ConstraintDiagnosticsReport`, + `summarise_residuals`, `UNKNOWN_EXPLANATION` or the misfit methods in + `docs/source` or `API.md` (checked 2026-11-05). + +#### S5 decisions + +1. **Resize contract.** For every structured grid, each assignment of `origin`, + `maximum`, `nsteps` or `step_vector` keeps the extent that the caller asked + for, or encloses it by less than one cell. `nsteps` counts nodes in the + geometry classes. A repeated assignment of the same value changes nothing. + Assignments of different values can add up. The sequence A, B, A can enlarge + the extent by less than one cell of B plus less than one cell of A. A test + measures this (`test_step_vector_a_b_a_extent_drift`). +2. **3D `nsteps` setter.** Keep origin and extent. The new step is + `extent / (nsteps_new - 1)`. The extent is `nsteps_cells * step_vector`, taken + before the change. A value below 2 on any axis raises `ValueError`. A value + that is not a whole number raises `ValueError`. The setter does not round, + because rounding hides a caller error. A float array of whole numbers is + accepted and stored as `int`. A scalar is broadcast, as before. +3. **`set_nelements` (Q4, answered 2026-11-05).** S5 accepts the exact extent + after `set_nelements`: `maximum` is the same before and after the call. The old + code changed it. `set_nelements` does `nsteps = np.maximum(np.ceil(...).astype(int), 2)`. + The clamp keeps behaviour for valid callers (a value of 2 or more is unchanged) + and avoids a new `ValueError` when the requested element count is very small. + The old code gave one node (zero cells) for such a request, which no caller + can use. The setter accepts only `nsteps >= 2`, because one node has no + extent to keep. +4. **3D `step_vector` setter.** Keep origin and the requested extent. The cell + count is `_cells_for_length(extent, step_vector=new_step)` and the node count + is that plus one. The stored step is a copy of the value that the caller gave. + The grid can end up to one cell larger than the extent, never smaller. The loop + in `set_interpolation_geometry` ends as before, because the cell count grows + when the step gets smaller. +5. **2D grid (Q1, answered 2026-11-05).** `nsteps` becomes a property. Its setter + coerces the value with `np.ceil(np.asarray(v)).astype(int)`, as `__init__` does + today. `nsteps_cells`, `n_cell_x` and `n_cell_y` become read-only properties + derived from `nsteps`. `maximum` becomes a read-only property with the legacy + value `origin + nsteps * step_vector`. This value counts nodes as cells. S5 does + not change it, because FR-028 asks to keep it. The four existing assertions + of `maximum` must pass with no change. After an `nsteps` assignment, the three + counts and `maximum` follow the new `nsteps` (before, they stayed stale). +6. **Rectilinear grid.** The `origin`, `nsteps` and `step_vector` setters raise + `NotImplementedError` after construction. Resize needs new node arrays, and + the grid cannot guess them. `maximum` already has no setter. +7. **Cell centres (FR-015).** No change to `BoundingBox.cell_centres` unless a new + test fails (OP12). S5 adds tests for rotated boxes and for the cell + association of the returned grid. If a test fails, the fix is the closed form + `origin + (arange(nsteps - 1) + 0.5) * node_step_vector` for each axis, in + place of `linspace`. The task records the failure first. +8. **Misfit methods (Q2, answered 2026-11-05).** The old public methods keep their + names and their public result shape. Two new methods return the corrected + result as a 1-D array: `evaluate_gradient_cosines` and `evaluate_value_residuals` + (both names are free; `git grep -n -w` on 2026-11-05 gives no match). The old + methods and the new methods call the same private code (one code path for each + quantity). That code works on copies, converts the stored local vectors with + `interpolator.local_to_world_vectors`, converts the stored local positions with + `bounding_box.reproject`, selects the columns `[:, :3]` (position) and + `[:, 3:6]` (vector), and evaluates through the feature. The in-place change, the + `grad[:, :3:6]` column error and the local-versus-world comparison are removed + from the old methods. +9. **Old result shape.** `evaluate_gradient_misfit` returns one row for each + non-empty family, as before: shape `(1, N)` for one family, and `(2, N)` for two + families of equal length. For two families of different length it returns one + row `(1, N1 + N2)`; before, this case raised `ValueError` (OP16). With no rows it + returns `np.array([])` of shape `(0,)`, as before. `evaluate_value_misfit` + returns shape `(N,)`, as before. +10. **Gradient value changes meaning.** The old value was the dot product of an + unnormalised model gradient and a column slice. The new value is the cosine + between the world data direction and the model gradient. This is a change of + behaviour of a public method. S5 records it in `COMPAT.md` as a migration + notice, and in `ROADMAP.md` as a decision entry (the constitution row for the + public API asks for it). The new name `evaluate_gradient_cosines` gives the + same numbers in one 1-D array. +11. **Stale model (Q3, answered 2026-11-05).** The misfit methods and + `verification_report` must not change state and must not rebuild. First they + ask `_stale_reason()`. It returns a text if `builder._up_to_date` is `False`, if + any fault builder of the builder is stale, or if `interpolator.up_to_date` is + `False`. It returns `None` otherwise (and when `builder` is `None`, only the + interpolator counts). If the text is not `None`, the method does not evaluate. + It writes one `warnings.warn(text, UserWarning, stacklevel=2)`. The text names + the stale object and says that the model must be built first. The new methods + then return `np.full(n, np.nan)`. The old methods return their usual shape + filled with NaN. Reasons: (a) the call changes no state (SC-009); (b) no + exception stops a loop over features; (c) NaN is the marker for an unusable row + in this design; (d) a number from old coefficients looks valid and is not; + (e) a `UserWarning` shows by default, a log line does not. When the model is up + to date, the methods evaluate with `self.evaluate_value` and + `self.evaluate_gradient`. Their call to `builder.up_to_date()` then does + nothing, because `_stale_reason()` has checked the same state (faults + included). The number `n` is the row count of the stored array. Reading a shape + changes nothing. +12. **Gradient misfit is a cosine.** The pure function `direction_cosines` returns + the cosine between the data direction and the model gradient for each row. A + row with a zero, not finite or missing vector gives `NaN`. +13. **Value misfit (Q5, answered 2026-11-05).** The public `evaluate_value_misfit` + keeps its old meaning: `|value - field|` divided by the field range + `self.max() - self.min()`. It divides only when the range is finite and greater + than zero. Otherwise it returns the unnormalised absolute difference. The + physical residual uses a private helper `_value_residuals_fresh()` that gives the + unnormalised `|value - field|` in field units. `evaluate_value_residuals`, + `evaluate_value_misfit` and the physical report entry `value` all call it. So + the unit `"field"` is true. +14. **Report.** `ConstraintDiagnosticsReport` gets three new fields with defaults, + appended after the existing ones. `to_dict` and `summary` show them only when + they are not empty. So a report that was built before S5 gives the same + output. A new method `GeologicalFeature.verification_report()` fills them. It + builds a new report. It does not store it in `latest_diagnostics_report`. Its + private call is `_build_constraint_diagnostics_report`. +15. **Separate residuals.** `physical_residuals` holds quantities in data units + (value: field units; direction: degrees). `assembled_residuals` holds the + `A x - b` rows of each constraint, in the assembler's normalisation. The two + dicts keep the plain names as keys. Assembled keys are constraint names (FDI: + `value`, `gradient`; PLI: `value`, `norm`). Physical keys are `value`, `gradient` + and `normal`. The two sets can share a key. S5 does not prefix the keys, + because the assembled keys must stay equal to the names that + `calculate_residual_for_constraints` uses. Each entry has a `kind` field + (`"physical"` or `"assembled"`). The two dicts are separate objects and never + merge in `to_dict`. The tests check the `kind` field and the object identity, + not a claim that the keys differ. +16. **Unsolved model.** `verification_report` fills the two residual blocks only if + `_stale_reason()` is `None` and `interpolator.c` is not `None` and the + interpolator has `calculate_residual_for_constraints`. If not, both blocks stay + empty (`{}`), one warning is written (stale case), and the exclusions are still + reported, because they read state only. `calculate_residual_for_constraints` + itself does not change in S5. +17. **Exclusions.** Each exclusion is an `ExclusionRecord` with a family, a reason, + a count and a plain explanation. Sources, in this order: + (a) The builder's `excluded_rows`. A builder record keeps its own constraint key + as `family` (for example `restricted_orthogonal`). It is not a data family. S5 + does not subtract it from any family count. + (b) The outside-model count of each family, reason `outside_model`, in points. + (c) For each family whose `dropped_rows` is not `None`: the rows that + `outside_model` does not explain, reason `unattributed`, in rows. The unit of + `dropped_rows` is the matrix row, so the outside count is first multiplied by + `rows_per_point` (2 for `gradient` when the interpolator is a + `DiscreteInterpolator`, 1 otherwise). The record exists when the difference is + positive. A family with `dropped_rows` equal to `None` (interface, + regularisation, other) gets no `unattributed` record, because it has no basis. + For Surfe, `dropped_rows = outside` and so the difference is zero. + A reason with no registered explanation keeps its count and gets + `UNKNOWN_EXPLANATION`. No row is dropped from the count. +18. **Out of S5.** Residual reports for fold cells and P2 cells (the cells that the + spec table links to FR-024 with SC-003, SC-007 and SC-008) are not in S5. S2 + covers the fold and P2 data paths. S6 covers the backend matrix. The summary + code does not depend on the engine, so S5 adds no cell for them. +19. **Python 3.9.** `_structured_grid_3d.py` uses `Optional[np.ndarray]` from + `typing`. It gets no `from __future__ import annotations` (the file stays + unchanged in style). `_diagnostics.py` already has the future import and uses + `float | None` in annotations. +20. **Exports.** `loop_interpolation.__init__` exports `ResidualSummary`, + `ExclusionRecord`, `direction_cosines` and `build_exclusion_records`. + `summarise_residuals`, `UNKNOWN_EXPLANATION` and `EXCLUSION_EXPLANATIONS` are + not exported. The feature code and the tests import them from + `loop_interpolation._diagnostics`. + +#### S5 contracts + +| Boundary | Input | Output | Owner | +|---|---|---|---| +| `StructuredGrid3DGeometry.nsteps` and `.step_vector` setters, `set_nelements` | node count, step, or element count | origin and extent kept (extent enclosed by less than one cell for the step setter) | `_structured_grid_3d.py` | +| `StructuredGrid2DGeometry.nsteps`, `maximum` and cell counts | `origin`, `nsteps`, `step_vector` | `nsteps` is an int array; counts derived, never stale; `maximum` legacy value `origin + nsteps * step_vector` | `_structured_grid_2d.py` | +| `RectilinearGrid` resize setters | any value | `NotImplementedError`, no state change | `_3d_rectilinear_grid.py` | +| `BoundingBox.cell_centres` | the box | local centres, one for each cell of `regular_grid(local=False)` (unchanged) | `_bounding_box.py` | +| `GeologicalFeature.evaluate_gradient_misfit` and `evaluate_value_misfit` | stored local constraint arrays | old public shape and names; world-correct; stored arrays, weights, coefficients and both up-to-date flags unchanged; stale model gives a warning and NaN | `_geological_feature.py` | +| `GeologicalFeature.evaluate_gradient_cosines` and `evaluate_value_residuals` | the same | 1-D arrays; same rules | `_geological_feature.py` | +| `direction_cosines` | two `(N, d)` arrays in one basis | `(N,)` cosines | `_diagnostics.py` | +| `GeologicalFeature.verification_report` | the feature | a new `ConstraintDiagnosticsReport` with residuals and exclusions; no state change | `_geological_feature.py` | + +#### S5 functions and signatures + +`StructuredGrid3DGeometry` (`_structured_grid_3d.py`; add `from typing import +Optional` at the top): + +```python +def _cells_for_length( + self, length: np.ndarray, step_vector: Optional[np.ndarray] = None +) -> np.ndarray: + """Number of cells per axis that cover ``length``. + + ``step_vector`` defaults to the current step. The rule of S1 is unchanged: + ``ceil(length / step - 1e-9)``; 0 cells become 2; a negative count raises + ``ValueError("Cannot resize the grid.")``. + """ +``` + +`set_nelements`: the line `nsteps = np.ceil(...).astype(int)` becomes + +```python +nsteps = np.maximum( + np.ceil((self.maximum - self.origin) / step_vector).astype(int), 2 +) +``` + +`nsteps` setter: + +```python +@nsteps.setter +def nsteps(self, nsteps): + nsteps = np.broadcast_to(np.asarray(nsteps), self._nsteps.shape) + if not np.all(np.isfinite(nsteps)) or np.any(nsteps != np.round(nsteps)): + raise ValueError("Cannot resize the grid. nsteps must be whole numbers.") + nsteps = nsteps.astype(int) + if np.any(nsteps < 2): + raise ValueError("Cannot resize the grid. nsteps must be at least 2.") + extent = self.nsteps_cells * self._step_vector + self._step_vector = extent / (nsteps - 1) + self._nsteps = np.array(nsteps) + self._on_geometry_change() +``` + +`step_vector` setter: + +```python +@step_vector.setter +def step_vector(self, step_vector): + step_vector = np.array(step_vector, dtype=float) + extent = self.nsteps_cells * self._step_vector + cells = self._cells_for_length(extent, step_vector=step_vector) + self._nsteps = cells + 1 + self._step_vector = step_vector + self._on_geometry_change() +``` + +`StructuredGrid2DGeometry` (`_structured_grid_2d.py`): in `__init__`, replace the +assignments of `self.nsteps`, `self.maximum`, `self.nsteps_cells`, `self.n_cell_x` +and `self.n_cell_y` (lines 23, 26, 28 to 30) with `self.nsteps = nsteps` (the setter +coerces) and keep the other lines. Add: + +```python +@property +def nsteps(self) -> np.ndarray: + return self._nsteps + +@nsteps.setter +def nsteps(self, nsteps): + self._nsteps = np.ceil(np.asarray(nsteps)).astype(int) + +@property +def nsteps_cells(self) -> np.ndarray: + return self.nsteps - 1 + +@property +def n_cell_x(self): + return self.nsteps[0] - 1 + +@property +def n_cell_y(self): + return self.nsteps[1] - 1 + +@property +def maximum(self) -> np.ndarray: + # Legacy value (FR-028): it counts nodes as cells, so it is one step + # beyond the last node. ``nodes`` uses ``nsteps_cells``. Do not change it. + return self.origin + self.nsteps * self.step_vector +``` + +The properties `n_cell_x` and `n_cell_y` return the same NumPy integer type as +the old attributes. The task greps for readers and writers of the five names +before the edit (`git grep -n "n_cell_x\|n_cell_y\|nsteps_cells\|\.maximum ="`). A +writer outside `__init__` must move to the derived form. The `__init__` order +matters: `nsteps` is set before `origin` is read. + +`RectilinearGrid` (`_3d_rectilinear_grid.py`): override the three setters. + +```python +_RESIZE_MESSAGE = ( + "A rectilinear grid has fixed node arrays and cannot be resized. " + "Build a new RectilinearGrid." +) + +@property +def origin(self): return self._origin +@origin.setter +def origin(self, origin): raise NotImplementedError(_RESIZE_MESSAGE) +# the same pattern for nsteps and step_vector +``` + +`_diagnostics.py` (new code; existing classes and fields keep their order): + +```python +@dataclass(frozen=True) +class ResidualSummary: + name: str # constraint or data family + kind: str # "physical" or "assembled" + unit: str # "field", "degree" or "assembled" + count: int # rows with a finite residual + nan_count: int # rows with no finite residual (not in the statistics) + rms: float | None # None when count == 0 + max_abs: float | None + mean_abs: float | None + +@dataclass(frozen=True) +class ExclusionRecord: + family: str # a data family, or the constraint key of a builder record + reason: str + count: int # points for outside_model, rows for unattributed + explanation: str + +EXCLUSION_EXPLANATIONS: dict[str, str] = { + "invalid_reference_value": "The reference value is NaN or not finite at the sampled point.", + "invalid_reference_direction": "The reference direction has zero length or is not finite.", + "evaluation_error": "The reference feature raised an error when it was evaluated.", + "outside_model": "The point lies outside the model bounding box. The count is in points.", + "unattributed": "The assembly has fewer rows than expected and no reason was recorded. The count is in matrix rows.", +} +UNKNOWN_EXPLANATION = "No explanation is registered for this reason." + +def summarise_residuals(name: str, kind: str, unit: str, residuals: np.ndarray) -> ResidualSummary: ... +def direction_cosines(model_vectors: np.ndarray, data_vectors: np.ndarray) -> np.ndarray: ... +def build_exclusion_records( + excluded_rows: dict[str, dict[str, int]], + families: dict[str, ConstraintFamilyDiagnostics], + rows_per_point: dict[str, int] | None = None, +) -> tuple[ExclusionRecord, ...]: ... +``` + +`summarise_residuals` flattens the array, counts non-finite values in `nan_count` +and computes `rms`, `max_abs`, `mean_abs` from the finite values. It does not +change its input. `direction_cosines` raises `ValueError` if the two shapes +differ. For each row it divides the dot product by the product of the two norms. +A row whose norm is zero or not finite gives `NaN`. Both inputs stay unchanged. + +`build_exclusion_records` returns the records sorted by `(family, reason)`. The +rules are: + +- One record for each reason in `excluded_rows`. The family is the key of + `excluded_rows`. The count is not subtracted from any family. +- One `outside_model` record for each family with `outside_model_point_count > 0`. +- One `unattributed` record for each family with `dropped_rows` not `None` and + `dropped_rows - outside_model_point_count * k > 0`. Here `k` is + `rows_per_point.get(family, 1)` (`rows_per_point` `None` means an empty dict). + The count of the record is the difference. +- A count of zero makes no record. + +The report class gets three fields after `outside_model_points`: + +```python +physical_residuals: dict[str, ResidualSummary] = field(default_factory=dict) +assembled_residuals: dict[str, ResidualSummary] = field(default_factory=dict) +exclusions: tuple[ExclusionRecord, ...] = () +``` + +`to_dict` adds the keys `"physical_residuals"`, `"assembled_residuals"` (a dict of +the fields of each summary) and `"exclusions"` (a list of dicts) only when the +field is not empty. `summary()` adds one line for each non-empty block, after the +existing lines. The properties `total_rows` and `active_families` do not change. + +`GeologicalFeature` (`_geological_feature.py`; the module has the future +import): + +```python +def evaluate_gradient_misfit(self) -> np.ndarray: + """Cosine between each stored direction and the model gradient (1 = equal). + + One row for each non-empty family (gradient, then normal): shape (1, N) for + one family, (2, N) for two of equal length, (1, N1 + N2) for two of + different length. Shape (0,) with no rows. The model gradient and the data + direction are compared in world coordinates. The method changes no stored + array, weight, coefficient or flag. A stale model gives a UserWarning and + NaN. The meaning of the value changed in S5 (see COMPAT.md). + """ + +def evaluate_value_misfit(self) -> np.ndarray: + """|value - field| at each stored value row, relative to the field range. + + The result is divided by the field range only when the range is finite and + greater than zero. Otherwise it is the absolute difference in field units. + """ + +def evaluate_gradient_cosines(self) -> np.ndarray: + """1-D cosines: gradient rows first, then normal rows. Shape (N,).""" + +def evaluate_value_residuals(self) -> np.ndarray: + """1-D |value - field| in field units, not normalised. Shape (N,).""" + +def verification_report(self) -> ConstraintDiagnosticsReport: + """Report with physical residuals, assembled residuals and exclusions.""" + +def _stale_reason(self) -> Optional[str]: ... +def _gradient_cosines_fresh(self) -> list: ... # one array for each non-empty family +def _value_residuals_fresh(self) -> np.ndarray: ... +def _physical_residual_summaries(self) -> dict[str, ResidualSummary]: ... +``` + +`_stale_reason` checks `self.builder._up_to_date`, then the same attribute of +every builder in `self.builder.faults`, then `self.interpolator.up_to_date`. It +returns a text for the first stale object, or `None`. + +Method body of `_gradient_cosines_fresh`, in order: (1) for each of the arrays +`get_gradient_constraints()` and `get_norm_constraints()` that has rows, `rows = +np.array(array, dtype=float, copy=True)`; (2) `positions = rows[:, :3]`, and +`reproject` them if `bounding_box` is not `None`; (3) `data = rows[:, 3:6]`, and +apply `self.interpolator.local_to_world_vectors(data)` when `bounding_box` is +not `None`; (4) `model = self.evaluate_gradient(positions)`; (5) append +`direction_cosines(model, data)`. The method returns the list (one array for each +non-empty family). `evaluate_gradient_cosines` calls `_stale_reason`; on a stale +model it warns and returns `np.full(n, np.nan)`; otherwise it returns +`np.concatenate` of the list (`np.zeros(0)` if the list is empty). The old +`evaluate_gradient_misfit` calls the same code and then builds the shape of +decision 9. The call `self.builder.up_to_date()` is removed from all of them. + +`_value_residuals_fresh`: copy `get_value_constraints()`, reproject the positions, +return `np.abs(rows[:, 3] - self.evaluate_value(positions))`. +`evaluate_value_misfit` returns the stale NaN array, or divides that result by +`self.max() - self.min()` when the range is finite and positive. + +`verification_report` body: `base = self.interpolator._build_constraint_diagnostics_report()` +(a new object; it does not set `latest_diagnostics_report`). If +`_stale_reason()` is not `None`, write the warning and skip the two residual +blocks. If `interpolator.c` is `None`, or the interpolator has no +`calculate_residual_for_constraints`, skip the assembled block. Otherwise +`assembled` = for each `(name, residual)` in +`calculate_residual_for_constraints().items()`, `summarise_residuals(name, +"assembled", "assembled", residual)`. Then `physical = +self._physical_residual_summaries()`; `excluded = getattr(self.builder, +"excluded_rows", {})`; `rows_per_point = {"gradient": 2}` if the interpolator is +a `DiscreteInterpolator`, else `None`. Return `dataclasses.replace(base, +physical_residuals=physical, assembled_residuals=assembled, +exclusions=build_exclusion_records(excluded, base.families, rows_per_point))`. +The physical summaries are `value` (from `_value_residuals_fresh`, unit +`"field"`), and `gradient` and `normal` (`degrees(arccos(clip(cosine, -1, 1)))` +of the list from `_gradient_cosines_fresh`, unit `"degree"`). So there is one +code path for each quantity. + +The names `ResidualSummary`, `ExclusionRecord`, `EXCLUSION_EXPLANATIONS`, +`UNKNOWN_EXPLANATION`, `summarise_residuals`, `direction_cosines`, +`build_exclusion_records`, `verification_report`, `evaluate_gradient_cosines`, +`evaluate_value_residuals`, `_stale_reason`, `_gradient_cosines_fresh`, +`_value_residuals_fresh` and `_physical_residual_summaries` are free (`git grep +-n -w` on 2026-11-05; `_build_constraint_diagnostics_report` exists and is used). + +#### S5 tests + +`packages/loop_common/tests/test_resize_paths.py` (new). Grids: `G1` origin +`(0,0,0)`, `nsteps (11,6,4)`, step `(1,2,3)`; `G2` origin `(100,-20,5)`, +`nsteps (5,5,5)`, step `(0.5,0.5,0.5)`; `G3` origin `(-1e4,2e3,-30)`, `nsteps +(21,3,8)`, step `(10,100,7.5)`. The 2D grids are `H1` origin `(0,0)`, `nsteps +(5,5)`, step `(1,1)`, and `H2` origin `(10,-5)`, `nsteps (8,4)`, step `(0.5,2)`. + +| Test function | Cells | Checks | +|---|---|---| +| `test_reapplying_same_value_10_times_changes_nothing_3d` | 3 grids x 4 setters (origin, maximum, nsteps, step_vector) = 12 | 10 repeats of each assignment leave `origin`, `maximum`, `step_vector` and `nsteps` equal to the start: counts exactly, lengths within `1e-12 * max(1, L)` (SC-011) | +| `test_nsteps_setter_keeps_extent` | 3 | `nsteps` doubled then halved: `maximum` equal to the start within `1e-12 * L`; node count equals the request | +| `test_step_vector_setter_encloses_extent` | 3 | step x 1.5: new `maximum` is at least the old one and less than one new step above it | +| `test_step_vector_a_b_a_extent_drift` | 3 | step A, then B = 1.5 A, then A again: the final `maximum` minus the start is at least 0 and less than `A + B` on each axis; record the measured value | +| `test_nsteps_below_two_is_rejected_without_change` | 1 | `ValueError`; state unchanged | +| `test_nsteps_not_whole_is_rejected_without_change` | 1 | `nsteps = (5.5, 5, 5)` and `NaN` raise `ValueError`; state unchanged; `(5.0, 5.0, 5.0)` is accepted and stored as `int` | +| `test_set_nelements_keeps_extent` | 3 | after `set_nelements(200)` `maximum` equals the start within `1e-12 * L`; the returned value equals `n_elements`; `nsteps` equals `max(ceil(extent / h), 2)` with `h = (volume / 200) ** (1 / 3)` computed in the test | +| `test_set_nelements_small_request_is_clamped` | 1 | `G1.set_nelements(1)` raises nothing; every `nsteps` is at least 2; the extent is kept | +| `test_2d_counts_are_derived` | 2 grids x 3 checks = 6 | after `nsteps` assignment, `nsteps_cells`, `n_cell_x`, `n_cell_y` and `n_elements` agree with the new `nsteps`; `nsteps` is an int array when the input is a float list; `maximum == origin + nsteps * step_vector` (legacy value, one step beyond `nodes.max(axis=0)`) | +| `test_rectilinear_resize_is_rejected` | 3 setters | `NotImplementedError`; nodes, origin and counts unchanged | + +`packages/loop_common/tests/test_cell_centre_sampling.py` (new). Boxes: `unit` +and `offset_anisotropic` from `test_bounding_box_structured_grid.py`. Frames: +rotations 0, 30 and 90 degrees about z; local origin `(0,0,0)` and `(7,-3,11)`. + +| Test function | Cells | Checks | +|---|---|---| +| `test_each_cell_has_exactly_one_centre` | 12 | `grid.cell_index_at(box.project(world_centres))` is `arange(n_cells)` with `grid = box.structured_grid(local_coordinates=True)`: no `-1`, no repeat, `np.bincount` all `1`. This is the main evidence for FR-015 and SC-011 | +| `test_local_coordinates_grid_associates_own_centres` | 12 | every `box.project(world_centres)` falls in a cell of `box.structured_grid(local_coordinates=True)`: `cell_index_at(...) != -1` for all points. The world centres are `reproject(box.cell_centres(order="F"))` | +| `test_centres_are_cell_midpoints_in_world` | 12 | `reproject(box.cell_centres(order="F"))` equals the mean of the 8 corner nodes of each cell of `box.structured_grid()`, within `Tp`. Not separate evidence for SC-011 (it restates the construction of `cell_centres`) | +| `test_centres_match_exported_cells` | 12 | world centres equal `grid.cell_centres` within `Tp`. Not separate evidence for SC-011 | +| `test_centres_inside_box` | 12 | each world centre is inside `[origin, maximum]` within `Tp`. Not separate evidence for SC-011 | + +The rotation is set with `box.set_local_transform(local_origin=..., local_rotation=...)`. +The task reads the S1 test `test_rigid_transform_contract.py` for the call and for +`Tp`. Rotation 0 with origin 0 is the control. Expected result by reading: all +pass on the current tree ("coverage, no defect", OP12). Only the first two tests +count as evidence for SC-011. + +`packages/loop_interpolation/tests/test_residual_exclusion_report.py` (new). It +tests the pure functions and the report class, with no model. + +| Test function | Cells | Checks | +|---|---|---| +| `test_direction_cosines_pure_function` | 6 | equal, opposite and orthogonal vectors; a zero row and a NaN row give `NaN`; shape mismatch raises `ValueError`; both inputs unchanged; scale invariance | +| `test_report_without_verification_is_unchanged` | 1 | `to_dict().keys()` equals the pre-S5 key set and the constructor still takes the old fields in order (rollback check) | +| `test_builder_records_keep_their_own_key` | 1 | `excluded_rows = {"restricted_orthogonal": {"invalid_reference_direction": 3}}` gives a record with `family == "restricted_orthogonal"` and count 3; no data-family count changes; no `unattributed` record appears | +| `test_unknown_reason_keeps_its_count` | 1 | a reason `"x"` in `excluded_rows` gives a record with `UNKNOWN_EXPLANATION` and the same count | +| `test_unattributed_rows_are_reported` | 1 | a `value` family with `dropped_rows` 5 and 2 outside points gives `unattributed` 3 | +| `test_gradient_outside_points_give_no_false_unattributed` | 1 | a `gradient` family with 10 source points, 2 outside points and `dropped_rows` 4 (`2 * 10 - 16` rows) and `rows_per_point={"gradient": 2}` gives no `unattributed` record; with `dropped_rows` 6 it gives 2 | +| `test_dropped_rows_none_gives_no_unattributed` | 1 | an `interface` family with `dropped_rows=None` and no outside points gives no record | +| `test_surfe_style_dropped_equals_outside` | 1 | `rows_per_point=None`, `dropped_rows` equal to the outside count: no `unattributed` record | + +`tests/unit/modelling/test_diagnostics_frames.py` (new). It uses +`world_foliation_model` (S4) and solves with the pattern of +`test_reference_constraint_frames.py`. + +| Test function | Cells | Checks | +|---|---|---| +| `test_misfit_leaves_state_unchanged` | 4 anchors x 3 rotations x 2 engines (FDI, PLI) = 24 | after a solve, snapshot (deep copies) of every `interpolator.data[...]` array, `c`, `interpolator.up_to_date`, `builder._up_to_date`, the weights of `constraints`, and `id(latest_diagnostics_report)`; call the two old misfit methods; every snapshot is equal (`np.array_equal`) (SC-009). It uses no new name | +| `test_new_misfit_methods_leave_state_unchanged` | 24 | the same snapshot; call `evaluate_gradient_cosines`, `evaluate_value_residuals` and `verification_report()` | +| `test_gradient_cosine_matches_analytic_direction_exact_field` | 24 | `monkeypatch` replaces `feature.evaluate_gradient` with a stub that returns the analytic world direction of the R1 foliation at the given points. R1 gradient and normal rows: `abs(1 - cos) <= 1e-10` for both the old method (one row) and `evaluate_gradient_cosines`. This tests the basis and column logic, not the solver (SC-009) | +| `test_gradient_cosine_matches_analytic_direction_solved_model` | 24 | the real solved model: `abs(1 - cos) <= tol_solved`. The value `tol_solved` is measured on the solved models. The task records the maximum and sets the bound with a margin of 2. It is not `1e-10` | +| `test_value_residuals_are_absolute_difference` | 4 anchors x 3 rotations = 12 | `evaluate_value_residuals()` equals `abs(value - evaluate_value(world_points))` at the 27 value rows; no NaN; unit field | +| `test_value_misfit_is_normalised_by_field_range` | 12 | `evaluate_value_misfit()` equals `evaluate_value_residuals() / (feature.max() - feature.min())` within `1e-12` | +| `test_value_misfit_with_zero_range_does_not_divide` | 2 | the model has value rows with one value (a near-constant field). The solver gives noise near `1e-12`, so the test forces the range with `monkeypatch`: `feature.max` and `feature.min` return the same number (cell 1) or `NaN` (cell 2). No `RuntimeWarning` (`warnings.simplefilter("error")`); the result equals the unnormalised absolute difference | +| `test_misfit_with_no_rows_is_empty` | 2 | no gradient rows: `shape == (0,)`; no value rows: `shape == (0,)` | +| `test_old_misfit_shapes_are_kept` | 4 | one family: `(1, N)`; two families of equal length: `(2, N)`; two of different length: `(1, N1 + N2)` (OP16); the values equal `evaluate_gradient_cosines()` | +| `test_stale_misfit_is_not_rebuilt` | 2 engines x 2 kinds (interpolator stale, builder stale) = 4 | uses the old methods only. Snapshot `interpolator.up_to_date`, `builder._up_to_date`, `c` and the data arrays. After a call, every snapshot is equal. On the current tree the call rebuilds the model, so `interpolator.up_to_date` changes from `False` to `True`: the test fails there (FR-031) | +| `test_stale_misfit_warns_once_and_returns_nan` | 4 | old and new methods each write exactly one `UserWarning` that names the stale object; the new methods return `np.full(n, np.nan)`; the old methods return the old shape filled with NaN | +| `test_stale_verification_report_does_not_rebuild` | 4 | `verification_report()` on a stale model: the state snapshot is equal; the two residual blocks are empty; the exclusions are present; one warning | +| `test_report_blocks_are_separate` | 2 engines = 2 | the two dicts are different objects; every physical entry has `kind == "physical"` and every assembled entry has `kind == "assembled"`; the physical `value` entry has unit `"field"` and its `rms` equals the RMS of `evaluate_value_residuals()` within `1e-12`; the assembled entries equal `summarise_residuals` of `calculate_residual_for_constraints()`; `to_dict` has two separate blocks | +| `test_unsolved_model_gives_no_residual_blocks` | 2 | a model with `interpolator.c = None`: the two blocks are empty and the exclusions are present; no exception | +| `test_exclusions_are_counted_and_explained` | 2 engines = 2 | the restricted-support reference case of S1 (`test_reference_constraint_frames.py`): the builder records appear under their constraint keys (`restricted_equality`, `restricted_orthogonal`) with the analytic count; the explanation is not empty and not `UNKNOWN_EXPLANATION` | +| `test_gradient_outside_points_give_no_false_unattributed_model` | 2 engines = 2 | a feature with gradient data and 3 points outside the model box: the report has an `outside_model` record of 3 and no `unattributed` record for `gradient` | + +Pre-fix evidence (FR-031). The state test with the old methods, the exact-field +cosine test and the stale-state test fail on the current tree. In the first, a +different array is stored in `data["gradient"]` after the call. In the second, +the cosine uses one column and local vectors. In the third, the call rebuilds the +model. The value-residual test records "not reproduced" if it passes on the +current tree (the value path does not change the stored positions). The resize +tests for `nsteps`, `step_vector`, `set_nelements` and the 2D counts fail on the +current tree (record the observed extent drift and the stale counts). The 2D +`maximum` is unchanged, so no 2D `maximum` test fails. The cell-centre tests are +coverage (OP12). The tests of the new names fail because the names do not exist. +The pre-fix part uses no new name. + +Existing tests that must pass with no change: +`packages/loop_common/tests/test_2d_discrete_support.py` and +`tests/unit/interpolator/test_2d_discrete_support.py` (the four `maximum` +assertions; S5 does not edit these two files), `test_rectilinear_grid.py`, +`test_bounding_box_structured_grid.py`, `tests/unit/modelling/test_scalar_field_grid.py`, +the constraint diagnostics tests (`git grep -ln "get_constraint_diagnostics_report" +tests packages`), and the S1 to S4 test files. + +Tests that change by design (Q4). The 3D golden test +`packages/loop_common/tests/test_structured_support_reference.py` and its data +file `packages/loop_common/tests/data/structured_support_reference.npz` change, +because the setters now keep the extent. The earlier sentence that this golden +test must pass with no change is withdrawn. The 2D golden test +`test_structured_grid_2d_reference.py` and `data/structured_grid_2d_reference.npz` +should not change (S5 keeps the 2D `maximum`). The task runs it, and if a key +differs, it records the key. The golden update follows this order: (1) run the +test and list each failing key; (2) review each key against the extent rule of +decision 1; (3) record each key, with the old and the new value, in the section +"S5 golden review (C-004)" of `verification.md` (SC-015: 0 tolerance changes); (4) +only then run `python packages/loop_common/tests/test_structured_support_reference.py +--update` (and the 2D script with `--update` only if a key differs); (5) run the +test again. + +Numeric-test survey (Q4). The tests that build a model run `set_nelements` and +the `step_vector` loop. The task runs `pytest tests/unit/modelling` and `pytest +packages/loop_interpolation/tests` on the fixed code, and `packages/loop_common/tests`. +For each test that fails or changes a number, it records the test, the old +value, the new value and the cause in "S5 golden review (C-004)". The directory +`packages/loop_interpolation/tests/golden` has its own option `--update-golden` +and is the parity oracle for the Rust path. If a file there changes, the task +records it and does not regenerate it before the review (OP17). + +#### S5 backend and capability cells + +| Cell | Result in S5 | +|---|---| +| FDI and PLI misfit, residual and exclusion report | Mandatory (the cells above). | +| P2 | Not in S5. The summary code does not depend on the engine. | +| Residual reports for fold and P2 cells (FR-024 with SC-003, SC-007, SC-008) | Not in S5. S2 covers the fold and P2 data paths, and S6 covers the matrix. | +| Native and Surfe | Not in S5. The misfit methods read the Python interpolator state only. Backend matrix: S6. | +| Resize of structured, 2D and rectilinear grids | Mandatory (resize tests above). Unstructured meshes cannot resize (OP4, S3). | + +#### S5 risks + +- The old misfit methods keep their public shape. The gradient value changes + meaning (decision 10). A caller outside the repository (the QGIS plugin, user + scripts) can depend on the old number. `COMPAT.md` and `ROADMAP.md` record it. + No caller in this repository uses the methods. +- `evaluate_gradient_misfit` used to call `builder.up_to_date()`, so it rebuilt a + stale model. Now it warns and returns NaN. A caller that relied on the rebuild + gets NaN and one warning, and must build the model first (OP13, closed). +- The extent change from `set_nelements` is a new risk. The old `nsteps` setter + kept `nsteps * step`. The new one keeps the extent. A model build runs + `set_nelements` and the `step_vector` loop, so grid steps, node coordinates and + solved values can change slightly in tests and in user models. The golden review + and the numeric-test survey record each change. The risk grows if an + interpolation golden file changes (OP17). +- Repeated `step_vector` assignments can add up to more than one cell of + enlargement (decision 1). A test measures the bound. +- The 2D `maximum` keeps the legacy value, so it still counts nodes as cells. A + reader that needs the true last node must use `nodes`. The value now follows + later `nsteps` assignments (it was stale before). +- `RectilinearGrid` now raises on resize. A caller that assigned `origin` or + `step_vector` got a wrong grid before. Now it gets a clear error. +- `evaluate_gradient` and `evaluate_value` of a feature call `builder.up_to_date()`. + The misfit methods call them only after `_stale_reason()` returned `None`. The + stale tests check this for FDI and PLI. +- The first S5 planner run read 14 code files, by targeted line ranges. The limit + of the spec-planner skill is 10 code files. The fix run read more. The design + uses only those reads. The risk is that an unread part of those files has other + readers. Each task greps for readers and writers before it edits. S5 changes 10 + code files at most (see "S5 files"). No task group changes more than 10. +- `build_exclusion_records` uses `rows_per_point` of 2 for the gradient family of + a `DiscreteInterpolator`. If an interpolator adds a different number of rows for + a gradient point, the `unattributed` count is wrong for that engine. The + gradient end-to-end test covers FDI and PLI only. + +#### S5 constitution check + +| Project rule | Result | How S5 agrees | +|---|---|---| +| Work stays within the existing packages | Pass | Three `loop_common` files, two `loop_interpolation` files (`_diagnostics.py` and `__init__.py`) and one `LoopStructural` file change. No new module. | +| Gates `pytest` and `ruff check LoopStructural packages` | Pass | The S5 row of "Stages and rollback" names both. | +| numpy docstrings, type hints, py39 | Pass | New functions and dataclasses have numpy docstrings and hints. `_diagnostics.py` has `from __future__ import annotations`. `_structured_grid_3d.py` uses `Optional`. | +| The public `GeologicalModel` API stays compatible; a break needs a ROADMAP.md decision and a major-version note | Pass | No `GeologicalModel` signature changes. The report only gets new defaulted fields. The misfit methods keep their names and public shape. The gradient value changes meaning: this is recorded in `COMPAT.md` (migration notice) and in a `ROADMAP.md` decision entry (Q2). The 2D `maximum` does not change. | +| Stages can be tested and reversed independently | Pass | S5 reverts to the S4 merge. | +| `ROADMAP.md` owns decisions; update it when one changes | Pass | S5 adds one decision entry for the misfit methods: "Misfit diagnostics of GeologicalFeature (spec 005, S5)". No other decision changes. | +| Feature documents stay under `specs/NNN-name/` | Pass | Only `plan.md`, `tasks.md` and `verification.md` change there. | + +#### S5 files + +At most 10 code files (.py) change in S5: 6 source files and 4 new test files. No +task group changes more than 10 code files. Data files and documents are not +counted. + +| File | Change | +|---|---| +| `packages/loop_common/src/loop_common/geometry/_structured_grid_3d.py` | `_cells_for_length` gets `step_vector`; new `nsteps` and `step_vector` setters; `set_nelements` clamp; `from typing import Optional` | +| `packages/loop_common/src/loop_common/geometry/_structured_grid_2d.py` | `nsteps` property with int coercion; `nsteps_cells`, `n_cell_x`, `n_cell_y` and `maximum` become properties (`maximum` keeps the legacy value, with a comment) | +| `packages/loop_common/src/loop_common/supports/_3d_rectilinear_grid.py` | `origin`, `nsteps`, `step_vector` setters raise `NotImplementedError` | +| `packages/loop_interpolation/src/loop_interpolation/_diagnostics.py` | `ResidualSummary`, `ExclusionRecord`, explanation table, three functions, three report fields, `to_dict` and `summary` | +| `packages/loop_interpolation/src/loop_interpolation/__init__.py` | Export four new names (`ResidualSummary`, `ExclusionRecord`, `direction_cosines`, `build_exclusion_records`). `summarise_residuals` and `UNKNOWN_EXPLANATION` are imported from `loop_interpolation._diagnostics` | +| `LoopStructural/modelling/features/_geological_feature.py` | Rewrite both misfit methods on shared private code; add `evaluate_gradient_cosines`, `evaluate_value_residuals`, `verification_report`, `_stale_reason` and the private helpers | +| `packages/loop_common/tests/test_resize_paths.py` | New | +| `packages/loop_common/tests/test_cell_centre_sampling.py` | New | +| `packages/loop_interpolation/tests/test_residual_exclusion_report.py` | New | +| `tests/unit/modelling/test_diagnostics_frames.py` | New | +| `packages/loop_common/tests/test_structured_support_reference.py` | Not edited. It must pass after the regeneration of its data file | +| `packages/loop_common/tests/data/structured_support_reference.npz` | Regenerated with `--update` after review (C-004, SC-015) | +| `packages/loop_common/tests/data/structured_grid_2d_reference.npz` | Regenerated only if the 2D reference test changes (not expected) | +| `packages/loop_common/tests/test_2d_discrete_support.py` and `tests/unit/interpolator/test_2d_discrete_support.py` | Not edited. They must pass with no change (Q1) | +| `COMPAT.md` | S5 notice: misfit gradient value (migration notice), stale model gives NaN and a warning, extent kept by the setters and `set_nelements`, rectilinear resize error, new report fields | +| `ROADMAP.md` | Decision "Misfit diagnostics of GeologicalFeature (spec 005, S5)" under "Decisions" | +| `specs/005-coordinate-frame-correctness/verification.md` | S5 pre-fix, post-fix, golden review (C-004), numeric-test survey, backend matrix, rollback and run-time sections | + +Not changed (the survey in "Later stages" listed them): `_bounding_box.py` +(`cell_centres`; conditional, OP12), `_structured_grid.py`, +`supports/_3d_structured_grid.py`, `_base_geological_feature.py`, +`_model_exporter.py`, `_discrete_interpolator.py` and `API.md`. `docs/source` and +`API.md` do not mention the changed names (`git grep`, 2026-11-05), so they need +no edit. + +The `COMPAT.md` notice says: `GeologicalFeature.evaluate_gradient_misfit()` +changed the stored constraint arrays in place, compared a local vector with a +world vector, and read one column. It now keeps its name and its public shape. It +changes nothing. It returns the cosine between the world data direction and the +model gradient for each row. The old value was a dot product with an unnormalised +model gradient. The new value differs. Use `evaluate_gradient_cosines()` for a +1-D array. The method no longer rebuilds a stale model: it writes a `UserWarning` +and returns NaN. `evaluate_value_misfit()` keeps its normalisation by the field +range. It no longer divides by zero. Assigning `nsteps` or `step_vector` on a 3D +grid, and `set_nelements`, keep the extent. A rectilinear grid refuses to resize. +`ConstraintDiagnosticsReport` has three new optional fields. With an identity +frame and no resize, the value misfit is the same as before. + +The `ROADMAP.md` entry gives the decision date `2026-11-05` (spec 005, S5, Q2). +It says that the misfit methods keep their names and public shape, that the +gradient value is now a cosine in world coordinates, that a stale model gives a +warning and NaN, and that the new names `evaluate_gradient_cosines` and +`evaluate_value_residuals` give the 1-D result. It lists the compatibility notice +in `COMPAT.md`. + +#### S5 task order (for the spec-tasks skill) + +No two tasks in one wave write `verification.md`. Two tasks never write the same +file in one wave. Each task group has at most 10 code files. + +1. The two `loop_common` test files (resize, cell centres), in two tasks that can + run in parallel. They run on the current tree. +2. The pre-fix evidence for step 1 in `verification.md` (one task, own wave). +3. The 3D setter fix with the `set_nelements` clamp, the 2D properties and the + rectilinear guard (three tasks, separate files, parallel). The 2D task greps + the readers and writers of the five 2D names first (`git grep -n + "n_cell_x\|n_cell_y\|nsteps_cells\|\.maximum ="`). +4. The numeric-test survey and golden review (one task, after step 3). It runs + `tests/unit/modelling`, `packages/loop_interpolation/tests` and + `packages/loop_common/tests`, records each change in "S5 golden review + (C-004)", and regenerates the golden `.npz` file after review. All later tasks + that build models depend on this step. It is the only `verification.md` writer + of its wave. +5. The post-fix evidence for the resize and cell-centre tests (one task, own wave, + after step 4). +6. `_diagnostics.py` (dataclasses, functions, report fields), with the pure-function + tests in `test_residual_exclusion_report.py`. Then the `__init__` export. This + step does not depend on steps 3 to 5 and can run in parallel with them, but its + evidence goes to `verification.md` in a later wave. +7. The feature tests in `test_diagnostics_frames.py` (state, cosine, value, empty, + shapes, stale, residual, exclusion), written before the misfit change. This task + owns that file alone. The pre-fix evidence goes to `verification.md` in its own + wave. It depends on step 4, because the solved-model tolerance and the values + depend on the built support. +8. The misfit rewrite and `verification_report` (one file; one task). It depends on + steps 6 and 7. +9. `COMPAT.md`, `ROADMAP.md`, the rollback check, and last the gates (the + orchestrator runs the full suites when memory is tight). `verification.md` is + written in the last wave only by this step. + +S5 rollback check. Revert to the S4 merge. Then run these checks. + +- (a) The field names and order of `ConstraintFamilyDiagnostics`, + `RegionCoverageDiagnostics` and the first five fields of + `ConstraintDiagnosticsReport` are equal before and after (`dataclasses.fields`). +- (b) For a report built with no verification, `to_dict()` is equal before and + after. +- (c) `from loop_interpolation import *` on the S4 code gives a subset of the S5 + names. +- (d) The signatures of `GeologicalFeature.evaluate_gradient_misfit` and + `evaluate_value_misfit` are equal. The result shape for one family is equal. +- (e) With an unchanged grid, the 3D `origin`, `maximum`, `nsteps` and + `step_vector` of a structured support are equal on both trees. +- (f) The four 2D `maximum` assertions pass on both trees. + +The result of each check goes to `verification.md`. + ### Stages and rollback | Stage | Files and behaviour | Acceptance check | Rollback target and compatibility check | @@ -1694,7 +2647,7 @@ are equal. | S2 | Fold geometry, directional regularisation, quadrature, masks, sign alignment for FDI-fold, PLI discrete fold, P2 and P2-fold; constant-norm variants (FR-008; R3 for P2; R4) ("Stage S2 design") | New S2 tests (SC-002 for constant norm; SC-007, SC-008 and SC-013 for every mandatory fold and P2 cell); the existing fold package tests pass with no change; `pytest`; `ruff check LoopStructural packages`; package tests | Revert to the S1 merge. Check: the S1 tests still pass. No saved format changed. | | S3 | Persistence format version, transform and coordinate-space metadata ("Stage S3 design"); no second projection at load; copies keep the frame; legacy load requires a transform or identity declaration; transform change policy and restored-point cache invalidation (FR-016 to FR-020); ROADMAP.md decision for FR-019; `COMPAT.md` migration note | SC-006, SC-012, SC-016, SC-017 | Revert to the S2 merge. Check: old-format R1 files load through their documented path; a new-format file fails with 1 actionable error and 0 state changes on the old reader (SC-017). | | S4 | Display, export of observation geometry, fault-centre geometry and intrusion classification use stored world locations (FR-022; "Stage S4 design") | World-location comparison within `Tp` (U3.4, U3.5) | Revert to the S3 merge. Check: viewer import paths unchanged. | -| S5 | Cell-centre sampling, remaining resize paths, non-mutating diagnostics, residual and exclusion report (FR-014 remainder, FR-015, FR-023, FR-024) | SC-009, SC-011 | Revert to the S4 merge. Check: diagnostics report fields unchanged for callers. | +| S5 | Remaining resize paths (3D `nsteps` and `step_vector`, `set_nelements`, 2D derived counts, rectilinear guard), golden review, cell-centre tests, non-mutating misfit diagnostics, residual and exclusion report (FR-014 remainder, FR-015, FR-023, FR-024; "Stage S5 design") | SC-009, SC-011; the new S5 tests; `pytest`; `ruff check LoopStructural packages`; package tests | Revert to the S4 merge. Check: diagnostics report fields unchanged for callers (`dataclasses.fields`, `to_dict` keys, exported names; "S5 rollback check"). | | S6 | Auxiliary point transformations and duplicate entry points (FR-025); displacement validation R2 (FR-026, FR-027); backend matrix incl. native and Surfe (FR-030); golden review (FR-029); 005 completion line in ROADMAP.md | SC-004, SC-010, SC-014, SC-015, SC-018 | Revert to the S5 merge. Check: `LoopStructural.utils` and `loop_common.math` import paths unchanged. | ## Test strategy @@ -1858,14 +2811,12 @@ grep survey on 2026-11-04 and are not final. `_3d_viewer.py`, `intrusion_builder.py`, `_fault_segment.py` and `exporters.py`; the code shows that S4 does not change them. S4 changes `_geological_feature.py` and `intrusion_frame_builder.py`. -- **S5**: `packages/loop_common/src/loop_common/geometry/_bounding_box.py` - (`cell_centres`), `_structured_grid.py`, `_structured_grid_2d.py`, - `packages/loop_common/src/loop_common/supports/_3d_structured_grid.py`, - `_3d_rectilinear_grid.py`, - `packages/loop_interpolation/src/loop_interpolation/_diagnostics.py`, and - `LoopStructural/modelling/features/_geological_feature.py` (found in the S4 - survey: the misfit methods change the stored arrays in place with - `grad /= norm`, and `grad[:, :3:6]` selects the position columns). +- **S5**: planned. See "Stage S5 design" and "S5 files". The survey listed + `_bounding_box.py`, `_structured_grid.py`, `supports/_3d_structured_grid.py` + and `_3d_rectilinear_grid.py`; the code shows that S5 changes only the + rectilinear grid among them (plus `_structured_grid_3d.py` and + `_structured_grid_2d.py`, which the survey did not list). `_diagnostics.py` and + `_geological_feature.py` are confirmed. - **S6**: `packages/loop_common/src/loop_common/math/_transformation.py`, `LoopStructural/utils/_transformation.py`, `packages/loop_interpolation/src/loop_interpolation/_native.py`, @@ -1889,8 +2840,8 @@ grep survey on 2026-11-04 and are not final. | FR-011 | Slip orthogonalisation; extent-direction row | SC-003, SC-013 | None | | FR-012 | Stored `_fault_buffer`, `_force_mesh_geometry`; factory override removed | SC-005; U1.5 test | None | | FR-013 | `local_enclosure`; `set_interpolation_geometry` | SC-005; enclosure unit test | None | -| FR-014 | `_cells_for_length`; setters | Resize-10-times test | Other supports S5, OP1 | -| FR-015 | S5 | S5 | OP1 | +| FR-014 | `_cells_for_length`; S1 setters; S5: 3D `nsteps` and `step_vector` setters, `set_nelements` clamp, 2D derived counts (legacy `maximum` kept), rectilinear guard | `test_resize_paths.py` (resize 10 times; extent kept; A, B, A drift; `set_nelements`); golden review (C-004) | Extent after `set_nelements` changes by design (Q4); OP17 | +| FR-015 | S5: tests for rotated boxes and cell association; `cell_centres` unchanged unless a test fails | `test_cell_centre_sampling.py` (two tests are evidence; three are coverage) | OP12 (defect not seen by reading) | | FR-016 | `frame_payload`; `to_dict` frame keys; format version 2 ("S3 saved format") | `test_frame_persistence.py::test_payload_declares_frame` | None | | FR-017 | `_restore_local_constraints` stores local arrays with no conversion | `test_local_arrays_unchanged_after_two_cycles` | None | | FR-018 | `_inherit_frame` in the three `copy()` methods | `test_frame_copy_and_lifecycle.py::test_copy_keeps_frame` | PLI has no `copy` (documented, not added); OP6 | @@ -1898,8 +2849,8 @@ grep survey on 2026-11-04 and are not final. | FR-020 | `frame_signature`, `_check_frame`, `StaleFrameError` (reject at use) | `test_stale_frame_rejected_without_state_change`, `test_fault_restoration_uses_current_frame` | Policy choice OP4; no restoration cache exists | | FR-021 | `validate_rigid_rotation` in `set_local_transform` and `rotation_xy`; saved metadata checked in `resolve_frame` | `test_rigid_transform_contract.py`; `test_malformed_metadata` | None | | FR-022 | Remove the extra `rescale` in `GeologicalFeature.get_data` and in the intrusion classification (`_network_model_values`); the fault centre is a world point since S1 and its consumers are unchanged | `test_world_location_consumers.py` (display, export, classification, 2D ellipse); S1 `test_fault_ellipsoid_centre_is_world` | OP8 to OP10 | -| FR-023 | S5 | S5 | OP1 | -| FR-024 | `excluded_rows` (S1 part) | Restricted-support exclusion count | Report integration S5, OP1 | +| FR-023 | S5: both misfit methods keep their names and shape, work on copies, compare in world coordinates, never rebuild (warn and NaN when stale); new `evaluate_gradient_cosines` and `evaluate_value_residuals`; `direction_cosines` | `test_diagnostics_frames.py`; `test_residual_exclusion_report.py` | OP16; the gradient value changes meaning (`COMPAT.md`, `ROADMAP.md`) | +| FR-024 | `excluded_rows` (S1); S5: `ResidualSummary`, `ExclusionRecord`, `build_exclusion_records`, report fields, `verification_report` | `test_exclusions_are_counted_and_explained`, `test_physical_and_assembled_residuals_are_separate` | Interpolator-side drops with no recorded reason are reported as `unattributed` (in rows). Residual reports for fold and P2 cells (SC-003, SC-007, SC-008) are out of S5: S2 and S6 | | FR-025 | S6 | S6 | OP1 | | FR-026 | S1: values, gradients, residuals (SC-003, SC-008); displacement S6 | `test_fault_frame_analytic.py` | Displacement OP1 | | FR-027 | S6 | S6 | OP1 | @@ -1916,9 +2867,9 @@ grep survey on 2026-11-04 and are not final. | SC-006 | Format version 2; local arrays restored with no conversion | `test_round_trip_world_values`, `test_local_arrays_unchanged_after_two_cycles` | None | | SC-007 | Reference sampling; S2 fold sampling and P2 | `test_reference_constraint_frames.py`; `test_fold_constraint_frames.py` | None | | SC-008 | Reference sampling; S2 fold rows and P2 | `test_reference_constraint_frames.py` (FDI, PLI, P2); `test_fold_constraint_frames.py` | None | -| SC-009 | S5 | S5 | OP1 | +| SC-009 | S5: misfit methods; `direction_cosines` | `test_misfit_leaves_state_unchanged`, `test_stale_misfit_is_not_rebuilt`, `test_gradient_cosine_matches_analytic_direction_exact_field`, `test_gradient_cosine_matches_analytic_direction_solved_model` (FDI and PLI; P2 not in S5) | None | | SC-010 | S6 | S6 | OP1 | -| SC-011 | `_cells_for_length` (resize part) | Resize test | Cell centres S5, OP1 | +| SC-011 | `_cells_for_length` and the S5 setters (resize part); S5 cell-centre tests | `test_resize_paths.py`, `test_cell_centre_sampling.py` (first two tests) | OP12, OP17 | | SC-012 | Rotation part in S1; metadata and legacy part in S3 | `test_rigid_transform_contract.py`; `test_legacy_*`, `test_malformed_metadata` | None | | SC-013 | Extent-direction row; S2 fold R4 directions | `test_fault_frame_analytic.py`; `test_fold_constraint_frames.py` | None | | SC-014 | Merge gates in every stage | `pytest`, ruff, package tests | None | @@ -1972,6 +2923,23 @@ grep survey on 2026-11-04 and are not final. display data through `get_data` and get world locations. - **Exclusion record**: `GeologicalFeatureBuilder.excluded_rows: dict[str, dict[str, int]]`. S5 connects it to the diagnostics report. +- **S5 diagnostics (owner: `loop_interpolation`)**: `ResidualSummary`, + `ExclusionRecord`, `direction_cosines(model_vectors, data_vectors)` and + `build_exclusion_records(excluded_rows, families)` are exported. + `ConstraintDiagnosticsReport` has the optional fields `physical_residuals`, + `assembled_residuals` and `exclusions`; `to_dict` shows them only when they + are not empty. `GeologicalFeature.verification_report()` fills them, makes + the private call `_build_constraint_diagnostics_report`, and does not store the + report. `summarise_residuals` and `UNKNOWN_EXPLANATION` are imported from + `loop_interpolation._diagnostics` and are not exported. The feature misfit + methods `evaluate_gradient_misfit` and `evaluate_value_misfit` keep their names + and public shape, work on copies and never rebuild. A stale model gives one + `UserWarning` and NaN. The gradient value is now a cosine in world coordinates. + The new methods `evaluate_gradient_cosines` and `evaluate_value_residuals` give + 1-D arrays. Spec 004 can report its directional rows through + `ResidualSummary(kind="assembled")`. Structured 3D grids keep origin and extent + on every setter and on `set_nelements`. The 2D `maximum` keeps its legacy value. + A rectilinear grid refuses to resize. - **Fault geometry**: an explicit `fault_center` is a world point; the generated centre extent row is `coord=2` with `strike_unit / major_axis`. Spec 003 golden fault models can change only under C-006/C-004 with a diff --git a/specs/005-coordinate-frame-correctness/run-state.json b/specs/005-coordinate-frame-correctness/run-state.json index 568a71c1..bcfee67b 100644 --- a/specs/005-coordinate-frame-correctness/run-state.json +++ b/specs/005-coordinate-frame-correctness/run-state.json @@ -1,10 +1,10 @@ { "spec": "005-coordinate-frame-correctness", - "stage": "S4", - "worktree": "C:/Users/lachl/dev/LoopStructural-005-s4", - "branch": "feat/005-s4-world-locations", + "stage": "S5", + "worktree": "C:/Users/lachl/dev/LoopStructural-005-s5", + "branch": "feat/005-s5-diagnostics-and-resize", "review_mode": "wave", - "start_commit": "9b3507a9a4d4c7d7a903a08c7d33619a7546bf61", + "start_commit": "1c029d21f92177e8d529d64b334c3cb0ca6ec29e", "runtime": { "runtime": "claude", "max_workers": 5, @@ -14,243 +14,250 @@ "review_spec": "opus", "review_verifier": "opus" }, - "env": "PYTHONPATH=C:/Users/lachl/dev/LoopStructural-005-s4;C:/Users/lachl/dev/LoopStructural-005-s4/packages/loop_common/src;C:/Users/lachl/dev/LoopStructural-005-s4/packages/loop_interpolation/src; python = main checkout .venv (C:/Users/lachl/dev/LoopStructural/.venv/Scripts/python.exe); run the root pytest and the package tests in separate calls, one heavy run at a time (low memory); detached script + foreground poll for long runs; logs with Out-File -Encoding utf8; for the viewer test (T054) add /packages/loopstructural_visualisation/src to PYTHONPATH (else it skips); use PowerShell for pytest (Bash resolves loop_interpolation from the main checkout); pyevtk for test_exporters: pip install --target /pyevtk_pkg pyevtk (system python), add the dir to PYTHONPATH", - "base_note": "Built from docs/005-s4-plan (9b3507a9) on origin/master (S1 to S3 and PR #400 merged). Not pushed; no PR.", + "env": "PYTHONPATH=C:/Users/lachl/dev/LoopStructural-005-s5;C:/Users/lachl/dev/LoopStructural-005-s5/packages/loop_common/src;C:/Users/lachl/dev/LoopStructural-005-s5/packages/loop_interpolation/src;C:/Users/lachl/dev/LoopStructural-005-s5/packages/loopstructural_visualisation/src;C:/Users/lachl/AppData/Local/Temp/claude/C--Users-lachl-dev-dev-tools/f206cbcd-e5d0-45be-a4c8-5f5c89bb26f1/scratchpad/pyevtk_pkg; python C:/Users/lachl/dev/LoopStructural/.venv/Scripts/python.exe; run pytest from PowerShell; root pytest and package tests in separate calls", + "base_note": "Built from docs/005-s5-plan (1c029d21), which sits on feat/005-s4-world-locations (S4, PR #401, not merged). Not pushed; no PR.", "waves": [ { - "wave": 23, - "base": "9b3507a9a4d4c7d7a903a08c7d33619a7546bf61", - "tasks": [ - "T051" - ], - "check": "scope pass; ruff pass; payload_declares 36 passed", - "commit": "7b4596046b230fb5f8cd5f55c40e6e277d85f932", - "review": "pass: sonnet, no findings" + "wave": 32, + "base": "e583f835337c5d23794f41e2fd6597c4064a0832", + "tasks": "T069,T070,T071", + "check": "scope pass; ruff pass; T069 24 failed/76 passed incl T070 (expected pre-fix); T071 ImportError collection (expected)", + "commit": "cc7d5f5ee0795bbec1d5aee6ece628841858fb4b", + "review": "fix round then recheck: pass (1 major fixed: T070 test premise/labels; 2 minor fixed); fix commit 04420f44" }, { - "wave": 24, - "base": "7b4596046b230fb5f8cd5f55c40e6e277d85f932", - "tasks": [ - "T052" - ], - "check": "scope pass; ruff pass; 18 failed 1 passed pre-fix (shift by local origin; rotation)", - "commit": "448aba2f29d1dfef4d1a9060335a6e94a010768d", - "review": "pass: sonnet, no findings" + "wave": 33, + "base": "436d5f3f9d282c6048c7b27c268acdd975f166fa", + "tasks": "T072,T073", + "check": "scope pass; residual report test 9 passed; ruff pass", + "commit": "0b4eef55d3e0896875495a3e783081b5a0cc6f84", + "review": "pass (no findings; 2 optional edge notes: finite**2 overflow near 1e154, inf*0 warning in direction_cosines)" }, { - "wave": 25, - "base": "448aba2f29d1dfef4d1a9060335a6e94a010768d", - "tasks": [ - "T053" - ], - "check": "scope pass; ruff pass; 12 failed AttributeError pre-fix", - "commit": "1b4abbb0020cce326ef56c7452de704af79513ad", - "review": "pass: sonnet, no findings" + "wave": 34, + "base": "6c7d80ecf8dabe9871156893783afbb8baacb00e", + "tasks": "T074,T075,T076,T077", + "check": "scope pass; resize+rectilinear+2d tests 78 passed; ruff pass", + "commit": "71a7948a4de89574ee28dded5b564e75970b4cf4", + "review": "pass (3 minor, open)" }, { - "wave": 26, - "base": "1b4abbb0020cce326ef56c7452de704af79513ad", - "tasks": [ - "T054" - ], - "check": "scope pass; ruff pass; fault_ellipse 4 passed (coverage); file collects without viewer", - "commit": "121b75c6a62c88188de06f6869cd132b35a712db", - "review": "pass: sonnet, 1 minor open" + "wave": 35, + "base": "b93c067182c78a0429e7934199eef890469e4bbf", + "tasks": "T078", + "check": "scope pass; golden tests 9 passed; survey: modelling 1540 passed, interpolation 615 passed, loop_common 406 passed, tests/unit rest 504 passed; 78 of 611 3D golden keys changed (causes C1 C2 C3); 2D golden unchanged; loop_interpolation golden unchanged", + "commit": "09d2e7b6bc82fcef0d813f45020015c171b12d23", + "review": "pass (no findings; reviewer confirmed 78 keys, causes C1 C2 C3)" }, { - "wave": 27, - "base": "27cd35bfc8692315a79ad6a06f019de6fc721ff8", - "tasks": [ - "T055", - "T056" - ], - "check": "scope pass; ruff pass; T056 25 failed (24 DID NOT WARN, 1 AttributeError); T052-T054 30 failed 5 passed", - "commit": "769562fe2ab0cb9caf32c3fc5c6653e5718d88cb", - "review": "pass: sonnet, 2 minor open" + "wave": 36, + "base": "8dff4ec4ce03031efb2e79e5a3b0322d1e8a7e92", + "tasks": "T079,T080", + "check": "scope pass; resize+cell-centre 100 passed; T080 127 failed pre-fix as expected; ruff pass", + "commit": "a69408bb24ad80e377edbc95472ff3a8f6c6a971", + "review": "pass (3 minor, open)" }, { - "wave": 28, - "base": "d8ed7d3eba533567046ca057963c7d98bf65b399", - "tasks": [ - "T057", - "T058", - "T059", - "T060", - "T061", - "T062" - ], - "check": "scope pass except approved T056 test fix; ruff pass; 62 passed (consumers + public api); test_exporters 12 passed 1 pre-existing fail (gocad)", - "commit": "7d5959e7f806ef8ff571406b0045be3c17301000", - "review": "pass: sonnet, no findings" + "wave": 37, + "base": "223c912bb05aab1ec23ca884fe8b19987ff53b3d", + "tasks": "T081", + "check": "scope pass; 20 new cells fail pre-fix as expected; ruff pass", + "commit": "88a7096f2d1c8f26e2abd86778cb08313d611ad4", + "review": "pass (3 minor, open; T084 context updated: guard both blocks on c is not None)" }, { - "wave": 29, - "base": "2e84bfcc443c6f60c9d2d55412f508ea0f0d227b", - "tasks": [ - "T063", - "T064", - "T065", - "T066" - ], - "check": "scope pass; greps ok; consumers 60 passed; intrusions 9; api contract 2; exporters 12 + 1 pre-existing fail", - "commit": "75c2cb3517b600036f847fbd7bf9af6e316ba583", - "review": "pass: sonnet, no findings" + "wave": 38, + "base": "e584bb20eb897241054ea2dda689211fdd48f1b8", + "tasks": "T082", + "check": "scope pass; 147 of 147 cells fail pre-fix; old method changes data[gradient] by up to 8.2254", + "commit": "6906c8ece48136d3dc16fede52876f01aa435195", + "review": "pass (no findings)" }, { - "wave": 30, - "base": "686322e9927b146c390947e562fe4388ab0549fa", - "tasks": [ - "T067" - ], - "check": "scope pass; (a)-(d) recorded; no rollback worktree left", - "commit": "0871a904d372076bb1e1ff5f62644598c44b8ba6", - "review": "pass: sonnet, no findings" + "wave": 39, + "base": "00c3e3741eb663e8660ebd44bc24fe987e37e2e1", + "tasks": "T083", + "check": "scope pass; 111 passed, 36 deselected; ruff pass", + "commit": "e97859bbb12fdc1d566cf5148e961b0874946c62", + "review": "pass (1 minor, open)" + }, + { + "wave": 40, + "base": "b8230d904fedce4ec53307afeeb788b039bbe43e", + "tasks": "T084", + "check": "scope pass; diagnostics+scalar_field 154 passed; ruff pass", + "commit": "0ef1d4396743856b803d6e0802b8dfa83504db83", + "review": "pass (no findings)" + }, + { + "wave": 41, + "base": "0764d785a5bac7e4253d067582309468a77a08f9", + "tasks": "T085,T086,T087", + "check": "scope pass; COMPAT greps ok; ROADMAP grep 1; rollback worktree removed; checks a to f pass", + "commit": "3a278aa5d439b60350fc5d350b4eb9a02cb297d5", + "review": "pass (no findings)" + }, + { + "wave": 42, + "base": "6b9e6467ff452970cb1c7fd3a39cc209cb9e2f36", + "tasks": "T088", + "check": "gates: ruff pass; packages 1021 passed 47 skipped 2 xfailed; root 2216 passed 4 skipped (936.67s)", + "commit": "pending", + "review": "covered by the end review" } ], "runs": { - "T051": [ + "T069": [ { "model": "sonnet", - "result": "pass: world_foliation_model added" + "result": "pass: 24 failing cells pre-fix" } ], - "T052": [ + "T070": [ { "model": "sonnet", - "result": "pass: 19 cells, 18 fail (shift by local origin), identity cell passes" + "result": "fail: test premise wrong (enclosure grid), 8 cells failed" + }, + { + "model": "sonnet", + "result": "pass: premise fixed, 60 cells pass; OP12 not reproduced (centres are correct)" } ], - "T053": [ + "T071": [ { "model": "sonnet", - "result": "pass: 12 cells fail with AttributeError" - }, + "result": "pass: ImportError as expected" + } + ], + "T072": [ { "model": "sonnet", - "result": "review fix pass: test_intrusion_steps_parameters_evaluate_model_at_world_points covers the production call site; mutation check fails as expected" + "result": "pass" } ], - "T054": [ + "T073": [ { "model": "sonnet", - "result": "pass: 4 cells coverage" + "result": "blocked: 4 T071 tests contradicted the task text" }, { "model": "sonnet", - "result": "fix pass: viewer import moved into the test with importorskip" + "result": "T071 tests corrected (fix worker); 9 passed" } ], - "T055": [ + "T074": [ { - "model": "sonnet", - "result": "pass: S4 pre-fix section in verification.md" + "model": "haiku", + "result": "pass" } ], - "T056": [ - { - "model": "sonnet", - "result": "pass: 25 failed (24 DID NOT WARN, 1 AttributeError _network_model_values)" - }, + "T075": [ { "model": "sonnet", - "result": "fix pass: scale inplace assertion records the existing behaviour (copy; caller array unchanged)" + "result": "pass: 29 passed" } ], - "T057": [ + "T076": [ { "model": "sonnet", - "result": "pass: 19 cells" + "result": "pass: 10+8+11 passed" } ], - "T058": [ + "T077": [ { "model": "haiku", - "result": "pass: 12 cells" + "result": "pass" } ], - "T059": [ + "T078": [ { - "model": "haiku", - "result": "fail: stacklevel=2 attributes the warning to the public_api wrapper; asked a question" - }, + "model": "sonnet", + "result": "pass" + } + ], + "T079": [ { - "model": "haiku", - "result": "pass after stacklevel=3 (fix); scale inplace test cell failed (test error, see answers)" + "model": "sonnet", + "result": "pass" } ], - "T060": [ + "T080": [ { - "model": "haiku", - "result": "pass: 4 rescale calls removed; 3 now-empty if real_coords blocks removed too; test_exporters ran only with pyevtk from a scratch --target dir" + "model": "sonnet", + "result": "pass: 127 failing cells pre-fix; tol_solved 9e-4 (measured max 4.414e-4)" } ], - "T061": [ + "T081": [ { - "model": "haiku", - "result": "pass" + "model": "sonnet", + "result": "pass: 20 failing cells pre-fix; PLI gradient unattributed 8 (k=2), outside_model 3, dropped_rows 14" } ], - "T062": [ + "T082": [ { - "model": "haiku", + "model": "sonnet", "result": "pass" } ], - "T063": [ + "T083": [ { "model": "sonnet", - "result": "pass" + "result": "pass: 111 passed" } ], - "T064": [ + "T084": [ { "model": "sonnet", - "result": "pass; Active written in the Removed-in column like the existing rows" + "result": "pass: 147 + 7 passed" } ], - "T065": [ + "T085": [ { - "model": "haiku", + "model": "sonnet", "result": "pass" } ], - "T066": [ + "T086": [ { "model": "haiku", "result": "pass" } ], - "T067": [ + "T087": [ { "model": "sonnet", - "result": "pass; Loop3DView.add_data does not exist: checked Loop2DView.add_data (signature equal) and Loop3DView.add_mesh" + "result": "pass: rollback commit 1c029d21 (S4 code); report has 4 fields on S4 not 5" } ], - "T068": [ + "T088": [ { "model": "sonnet", - "result": "pass; gates run by the orchestrator (detached script); Done-when root pytest has 1 pre-existing failure (gocad, not S4)" + "result": "pass (gates run by orchestrator)" } ] }, "findings": { - "found": 7, + "found": 19, "rejected": 0, "dropped": 0, - "fixed": 1, + "fixed": 4, "open_minor": [ - "T054 test_2d_viewer_fault_ellipse_centre_is_world: compares the patch centre to fault.fault_centre[:2] only; add a check that fault_centre[:2] is close to anchor[:2] (optional).", - "T056 test_internal_callers_do_not_warn: importorskip(pyevtk) sits after the first checks; move it to the start or split the exporter part.", - "T056 failing ids before the fixes (25: 24 DID NOT WARN, 1 AttributeError _network_model_values) are not in verification.md; record them in T063.", - "T067/T064: the plan names Loop3DView.add_data, which does not exist (Loop2DView.add_data does). COMPAT.md S4 notice repeats the name; correct in a later pass.", - "T054 fault ellipse test: add check against anchor[:2] (end review, minor).", - "T056 test_internal_callers_do_not_warn: importorskip(pyevtk) mid-test hides the checks that ran; move the exporter part to its own test (end review, minor).", - "T068 verification.md S4 deprecated-call grep: line numbers of verification.md hits are stale (wave 31, minor)." + "T075 _structured_grid_3d.py has a UTF-8 BOM added at line 1 (remove the 3 bytes).", + "T076 _structured_grid_2d.py has format-only hunks in is_inside, neighbour_global_indexes, cell_corner_indexes (noise).", + "T075 set_nelements passes a cell count to the nsteps setter, which takes nodes: nelements=1000 on a unit cube gives 9^3=729 cells (task text and old behaviour; not a regression).", + "T079 section of verification.md does not say that test_each_cell_has_exactly_one_centre was replaced, and calls the midpoint test evidence although it is a guard (OP12 not reproduced). T088 adds one sentence.", + "T080 exact-field stub ignores its points argument, so it cannot catch a wrong position frame (other tests catch it).", + "T080 one TOL_SOLVED=9e-4 for both engines (FDI max 1.2e-4 is looser).", + "T081 test_gradient_outside_points_with_pli_are_measured prints but does not assert unattributed == 8.", + "T081 _snapshot compares id(latest_diagnostics_report): hold the object and compare with is.", + "T083 _stale_reason checks one level of faults; BaseBuilder.up_to_date recurses, so a stale fault nested in a fresh fault builder is still rebuilt on evaluation (rare; plan names one level).", + "Spec review: _structured_grid_3d.py BOM (same as wave 34 minor).", + "Spec review: COMPAT.md ConstraintDiagnosticsReport row has the sentence about value misfit with identity frame and no resize; move it to the evaluate_value_misfit row (a built model is always resized).", + "Spec review: verification.md S5 test run time (OP2) count sentence: 150 more passed tests vs 147 cells + 1 S4 failure.", + "MAJOR OPEN (needs maintainer decision): _value_residuals_fresh and _gradient_cosines_fresh evaluate through the feature (plan decision 8), which applies faults a second time to stored positions that were already fault-restored. Physical residuals of a faulted feature compare two frames (measured: 0.080 vs true 0.0004 to 0.0018). Proposed fix: evaluate the interpolator directly at the stored restored positions with no fault step." ] }, - "answers": [ - "2026-11-05 (wave 28, orchestrator): (1) T059 stacklevel is 3, not 2, because @public_api adds one frame; the test requires the warning to name the caller file. (2) GeologicalModel.scale has always copied its input (np.array(points).astype(float)), so inplace=True never changed the caller array; FR-033 says keep behaviour, so the T056 scale test now checks the existing behaviour (result equals project(points); caller array unchanged) instead of mutation. (3) pyevtk is not in the venv: tests/unit/io/test_exporters.py and the exporter part of test_internal_callers_do_not_warn skip without it. pyevtk 1.7.0 was installed with pip --target into the scratch dir and put on PYTHONPATH for the checks; the venv is unchanged. (4) test_write_vol_gocad_is_broken_for_real_bounding_box fails before and after S4 (pre-existing: BoundingBox.__getitem__ raises LoopValueError for a slice; the test expects pytest.raises(LoopValueError) but the failure is at the base too) - to check." - ], + "answers": "2026-11-05 (S5 start): five gaps settled in plan.md Answers (bounding box = interpolator.bounding_box; stale recipe = assign up_to_date/_up_to_date False; 2D n_elements exists; report key set is a literal; PLI add_gradient_constraints is pass so PLI gradient family reports unattributed = 2*source_points).", "history": { "S1": { "stage": "S1", @@ -1179,19 +1186,303 @@ "review": "first review: 1 major (warning location and timing) + 2 minor; fixed with warn_legacy_identity() and 7 new test cells (commit 169ae619); fix checked by the orchestrator, no second reviewer round" } ] + }, + "S4": { + "spec": "005-coordinate-frame-correctness", + "stage": "S4", + "worktree": "C:/Users/lachl/dev/LoopStructural-005-s4", + "branch": "feat/005-s4-world-locations", + "review_mode": "wave", + "start_commit": "9b3507a9a4d4c7d7a903a08c7d33619a7546bf61", + "runtime": { + "runtime": "claude", + "max_workers": 5, + "worker_S": "haiku", + "worker_M": "sonnet", + "review_wave": "sonnet", + "review_spec": "opus", + "review_verifier": "opus" + }, + "env": "PYTHONPATH=C:/Users/lachl/dev/LoopStructural-005-s4;C:/Users/lachl/dev/LoopStructural-005-s4/packages/loop_common/src;C:/Users/lachl/dev/LoopStructural-005-s4/packages/loop_interpolation/src; python = main checkout .venv (C:/Users/lachl/dev/LoopStructural/.venv/Scripts/python.exe); run the root pytest and the package tests in separate calls, one heavy run at a time (low memory); detached script + foreground poll for long runs; logs with Out-File -Encoding utf8; for the viewer test (T054) add /packages/loopstructural_visualisation/src to PYTHONPATH (else it skips); use PowerShell for pytest (Bash resolves loop_interpolation from the main checkout); pyevtk for test_exporters: pip install --target /pyevtk_pkg pyevtk (system python), add the dir to PYTHONPATH", + "base_note": "Built from docs/005-s4-plan (9b3507a9) on origin/master (S1 to S3 and PR #400 merged). Not pushed; no PR.", + "waves": [ + { + "wave": 23, + "base": "9b3507a9a4d4c7d7a903a08c7d33619a7546bf61", + "tasks": [ + "T051" + ], + "check": "scope pass; ruff pass; payload_declares 36 passed", + "commit": "7b4596046b230fb5f8cd5f55c40e6e277d85f932", + "review": "pass: sonnet, no findings" + }, + { + "wave": 24, + "base": "7b4596046b230fb5f8cd5f55c40e6e277d85f932", + "tasks": [ + "T052" + ], + "check": "scope pass; ruff pass; 18 failed 1 passed pre-fix (shift by local origin; rotation)", + "commit": "448aba2f29d1dfef4d1a9060335a6e94a010768d", + "review": "pass: sonnet, no findings" + }, + { + "wave": 25, + "base": "448aba2f29d1dfef4d1a9060335a6e94a010768d", + "tasks": [ + "T053" + ], + "check": "scope pass; ruff pass; 12 failed AttributeError pre-fix", + "commit": "1b4abbb0020cce326ef56c7452de704af79513ad", + "review": "pass: sonnet, no findings" + }, + { + "wave": 26, + "base": "1b4abbb0020cce326ef56c7452de704af79513ad", + "tasks": [ + "T054" + ], + "check": "scope pass; ruff pass; fault_ellipse 4 passed (coverage); file collects without viewer", + "commit": "121b75c6a62c88188de06f6869cd132b35a712db", + "review": "pass: sonnet, 1 minor open" + }, + { + "wave": 27, + "base": "27cd35bfc8692315a79ad6a06f019de6fc721ff8", + "tasks": [ + "T055", + "T056" + ], + "check": "scope pass; ruff pass; T056 25 failed (24 DID NOT WARN, 1 AttributeError); T052-T054 30 failed 5 passed", + "commit": "769562fe2ab0cb9caf32c3fc5c6653e5718d88cb", + "review": "pass: sonnet, 2 minor open" + }, + { + "wave": 28, + "base": "d8ed7d3eba533567046ca057963c7d98bf65b399", + "tasks": [ + "T057", + "T058", + "T059", + "T060", + "T061", + "T062" + ], + "check": "scope pass except approved T056 test fix; ruff pass; 62 passed (consumers + public api); test_exporters 12 passed 1 pre-existing fail (gocad)", + "commit": "7d5959e7f806ef8ff571406b0045be3c17301000", + "review": "pass: sonnet, no findings" + }, + { + "wave": 29, + "base": "2e84bfcc443c6f60c9d2d55412f508ea0f0d227b", + "tasks": [ + "T063", + "T064", + "T065", + "T066" + ], + "check": "scope pass; greps ok; consumers 60 passed; intrusions 9; api contract 2; exporters 12 + 1 pre-existing fail", + "commit": "75c2cb3517b600036f847fbd7bf9af6e316ba583", + "review": "pass: sonnet, no findings" + }, + { + "wave": 30, + "base": "686322e9927b146c390947e562fe4388ab0549fa", + "tasks": [ + "T067" + ], + "check": "scope pass; (a)-(d) recorded; no rollback worktree left", + "commit": "0871a904d372076bb1e1ff5f62644598c44b8ba6", + "review": "pass: sonnet, no findings" + } + ], + "runs": { + "T051": [ + { + "model": "sonnet", + "result": "pass: world_foliation_model added" + } + ], + "T052": [ + { + "model": "sonnet", + "result": "pass: 19 cells, 18 fail (shift by local origin), identity cell passes" + } + ], + "T053": [ + { + "model": "sonnet", + "result": "pass: 12 cells fail with AttributeError" + }, + { + "model": "sonnet", + "result": "review fix pass: test_intrusion_steps_parameters_evaluate_model_at_world_points covers the production call site; mutation check fails as expected" + } + ], + "T054": [ + { + "model": "sonnet", + "result": "pass: 4 cells coverage" + }, + { + "model": "sonnet", + "result": "fix pass: viewer import moved into the test with importorskip" + } + ], + "T055": [ + { + "model": "sonnet", + "result": "pass: S4 pre-fix section in verification.md" + } + ], + "T056": [ + { + "model": "sonnet", + "result": "pass: 25 failed (24 DID NOT WARN, 1 AttributeError _network_model_values)" + }, + { + "model": "sonnet", + "result": "fix pass: scale inplace assertion records the existing behaviour (copy; caller array unchanged)" + } + ], + "T057": [ + { + "model": "sonnet", + "result": "pass: 19 cells" + } + ], + "T058": [ + { + "model": "haiku", + "result": "pass: 12 cells" + } + ], + "T059": [ + { + "model": "haiku", + "result": "fail: stacklevel=2 attributes the warning to the public_api wrapper; asked a question" + }, + { + "model": "haiku", + "result": "pass after stacklevel=3 (fix); scale inplace test cell failed (test error, see answers)" + } + ], + "T060": [ + { + "model": "haiku", + "result": "pass: 4 rescale calls removed; 3 now-empty if real_coords blocks removed too; test_exporters ran only with pyevtk from a scratch --target dir" + } + ], + "T061": [ + { + "model": "haiku", + "result": "pass" + } + ], + "T062": [ + { + "model": "haiku", + "result": "pass" + } + ], + "T063": [ + { + "model": "sonnet", + "result": "pass" + } + ], + "T064": [ + { + "model": "sonnet", + "result": "pass; Active written in the Removed-in column like the existing rows" + } + ], + "T065": [ + { + "model": "haiku", + "result": "pass" + } + ], + "T066": [ + { + "model": "haiku", + "result": "pass" + } + ], + "T067": [ + { + "model": "sonnet", + "result": "pass; Loop3DView.add_data does not exist: checked Loop2DView.add_data (signature equal) and Loop3DView.add_mesh" + } + ], + "T068": [ + { + "model": "sonnet", + "result": "pass; gates run by the orchestrator (detached script); Done-when root pytest has 1 pre-existing failure (gocad, not S4)" + } + ] + }, + "findings": { + "found": 7, + "rejected": 0, + "dropped": 0, + "fixed": 1, + "open_minor": [ + "T054 test_2d_viewer_fault_ellipse_centre_is_world: compares the patch centre to fault.fault_centre[:2] only; add a check that fault_centre[:2] is close to anchor[:2] (optional).", + "T056 test_internal_callers_do_not_warn: importorskip(pyevtk) sits after the first checks; move it to the start or split the exporter part.", + "T056 failing ids before the fixes (25: 24 DID NOT WARN, 1 AttributeError _network_model_values) are not in verification.md; record them in T063.", + "T067/T064: the plan names Loop3DView.add_data, which does not exist (Loop2DView.add_data does). COMPAT.md S4 notice repeats the name; correct in a later pass.", + "T054 fault ellipse test: add check against anchor[:2] (end review, minor).", + "T056 test_internal_callers_do_not_warn: importorskip(pyevtk) mid-test hides the checks that ran; move the exporter part to its own test (end review, minor).", + "T068 verification.md S4 deprecated-call grep: line numbers of verification.md hits are stale (wave 31, minor)." + ] + }, + "answers": [ + "2026-11-05 (wave 28, orchestrator): (1) T059 stacklevel is 3, not 2, because @public_api adds one frame; the test requires the warning to name the caller file. (2) GeologicalModel.scale has always copied its input (np.array(points).astype(float)), so inplace=True never changed the caller array; FR-033 says keep behaviour, so the T056 scale test now checks the existing behaviour (result equals project(points); caller array unchanged) instead of mutation. (3) pyevtk is not in the venv: tests/unit/io/test_exporters.py and the exporter part of test_internal_callers_do_not_warn skip without it. pyevtk 1.7.0 was installed with pip --target into the scratch dir and put on PYTHONPATH for the checks; the venv is unchanged. (4) test_write_vol_gocad_is_broken_for_real_bounding_box fails before and after S4 (pre-existing: BoundingBox.__getitem__ raises LoopValueError for a slice; the test expects pytest.raises(LoopValueError) but the failure is at the base too) - to check." + ], + "final": { + "head_gated": "1ba71f17", + "lint": "ruff check LoopStructural packages: pass", + "test_packages": "912 passed, 47 skipped, 2 xfailed in 150.24 s", + "test_root": "1 failed, 2066 passed, 4 skipped in 805.50 s; the failure test_write_vol_gocad_is_broken_for_real_bounding_box is pre-existing (confirmed with the S4 changes stashed; two LoopValueError classes); it was skipped before because pyevtk is missing from the venv (pyevtk from a scratch dir on PYTHONPATH)", + "format": "black --check: 3 library files dirty before S4 (exporters.py, intrusion_frame_builder.py, geological_model.py); only task hunks changed" + }, + "spec_review": { + "range": "9b3507a9..e0b52cc5", + "reviewer": "opus", + "verdict": "changes needed: 1 major (call site of _network_model_values not covered), 2 minor", + "fix": "1 major fixed with a new test (mutation proved)" + } } }, - "final": { - "head_gated": "1ba71f17", - "lint": "ruff check LoopStructural packages: pass", - "test_packages": "912 passed, 47 skipped, 2 xfailed in 150.24 s", - "test_root": "1 failed, 2066 passed, 4 skipped in 805.50 s; the failure test_write_vol_gocad_is_broken_for_real_bounding_box is pre-existing (confirmed with the S4 changes stashed; two LoopValueError classes); it was skipped before because pyevtk is missing from the venv (pyevtk from a scratch dir on PYTHONPATH)", - "format": "black --check: 3 library files dirty before S4 (exporters.py, intrusion_frame_builder.py, geological_model.py); only task hunks changed" + "s5_notes": { + "T070": "The task text premise (one centre per cell of structured_grid(local_coordinates=True)) is wrong for rotated boxes: that grid is the axis-aligned enclosure with a different step. Test rewritten on the box own cell grid. _bounding_box.py not changed. OP12: not reproduced.", + "T070_arg": "set_local_transform takes rotation_matrix=, not local_rotation=." + }, + "s5_notes2": { + "T071_T073": "Four T071 tests asserted () or no record where the T073 rule (outside_model record per family with outside points; builder keys are not families) gives outside_model records or needs family dropped_rows 0. Tests corrected to assert no unattributed record; implementation follows the task text." + }, + "s5_notes3": { + "T076": "Builder lines _geological_feature_builder.py 520/544/604 assign support.maximum; 3D supports only (2D supports are not built by the feature builder). 2D maximum is now read-only.", + "T075": "Writers checked: none depends on the old product rule; set_interpolation_geometry loop still terminates." + }, + "s5_notes4": { + "T080": "Test new_misfit_methods_leaves_state_unchanged name spelled with leaves so the -k filter matches. Contract: evaluate_gradient_cosines is flat (gradient rows then normal rows); old misfit ravel equals it; new method must call self.evaluate_gradient. TOL_SOLVED=9e-4." + }, + "s5_notes5": { + "T081": "PLI measured: 7 gradient rows stored (4 inside, 3 outside), builder dropped_rows=14, outside_model=3, unattributed=8 with rows_per_point gradient 2. FDI dropped_rows=6; needs rows_per_point={gradient:2} (else unattributed 3). Warning text must contain the word interpolator or builder." }, "spec_review": { - "range": "9b3507a9..e0b52cc5", + "range": "1c029d21..8ecb66fc", "reviewer": "opus", - "verdict": "changes needed: 1 major (call site of _network_model_values not covered), 2 minor", - "fix": "1 major fixed with a new test (mutation proved)" + "verdict": "changes needed: 1 blocker (confirmed by opus verifier: _stale_reason treats every fault as stale), 1 major (faulted feature compares frames twice; contradicts plan decision 8, NOT fixed, needs maintainer decision), 3 minor", + "fix": "blocker fixed by sonnet worker plus 2 regression tests; worker also fixed a loop-variable overwrite in evaluate_gradient of faulted features" + }, + "final": { + "head_gated": "8eb22671", + "lint": "ruff check LoopStructural packages: pass", + "test_packages": "1021 passed, 47 skipped, 2 xfailed in 98.44 s", + "test_root": "2218 passed, 4 skipped in 864.77 s", + "spec_review_recheck": "pass (opus), no findings" } } \ No newline at end of file diff --git a/specs/005-coordinate-frame-correctness/tasks.md b/specs/005-coordinate-frame-correctness/tasks.md index 7b85698e..713e2877 100644 --- a/specs/005-coordinate-frame-correctness/tasks.md +++ b/specs/005-coordinate-frame-correctness/tasks.md @@ -1,4 +1,4 @@ -# Tasks: Coordinate and fault-frame correctness (stages S1 to S4) +# Tasks: Coordinate and fault-frame correctness (stages S1 to S5) **Feature**: `005-coordinate-frame-correctness` | **Plan**: [plan.md](plan.md) @@ -10,15 +10,15 @@ Rules: - [P] = can be done at the same time as the other [P] tasks of the same wave. - Written in ASD-STE100 Simplified Technical English. - This file covers stages S1 (T001 to T017, done), S2 (T018 to T032), S3 - (T033 to T050, done) and S4 (T051 to T068) of plan.md. Stages S5 and S6 wait - for OP1 (a later spec-planner run for each stage, then a new spec-tasks run). + (T033 to T050, done) and S4 (T051 to T068, done) and S5 (T069 to T088) of plan.md. Stage S6 waits + for OP1 (a later spec-planner run, then a new spec-tasks run). --> ## Orchestrator notes - Run this file with the spec-orchestrator skill. It has the rules for the worktree, the sub agents, the checks and the commits. -- Waves: 31 (S1: 1 to 7, done; S2: 8 to 13; S3: 14 to 22, done; S4: 23 to 31). Tiers: [S] = simple model, [M] = standard model. +- Waves: 42 (S1: 1 to 7, done; S2: 8 to 13; S3: 14 to 22, done; S4: 23 to 31, done; S5: 32 to 42). Tiers: [S] = simple model, [M] = standard model. - Do not start a task with "Blocked by". ## Phase 1: Setup @@ -1166,6 +1166,221 @@ frame. A legacy file needs an explicit choice. A frame change never gives a stal Add the S4 gate counts to the post-fix section. Add to the status log of `ROADMAP.md` a line dated today: "Spec 005 stage S4 implemented: display, export and intrusion classification use stored world locations; rescale and scale deprecated." Report any library warning with the text "is deprecated" in the pytest summary. Write in ASD-STE100 Simplified Technical English. +## Phase 18: Stage S5 tests (FR-014, FR-015, FR-023, FR-024) + +**Goal**: Tests for the remaining resize paths, the cell-centre sampling and the residual and exclusion report are written before the code that they test. +**Independent test**: `pytest packages/loop_common/tests/test_resize_paths.py packages/loop_common/tests/test_cell_centre_sampling.py` and `pytest packages/loop_interpolation/tests/test_residual_exclusion_report.py`. + +- [x] T069 [US4] [M] Test: the resize paths keep the extent and the counts (U4.2; FR-014; SC-011). + - Files: `packages/loop_common/tests/test_resize_paths.py` + - Needs: T068 + - Wave: 32 + - Done when: `pytest packages/loop_common/tests/test_resize_paths.py` fails before T075, T076 and T077 (record the failing ids and the first error of each) and passes after T077. + - Context: New file. Classes: `StructuredGrid3DGeometry` (`loop_common.geometry._structured_grid_3d`), `StructuredGrid2DGeometry` (`_structured_grid_2d`), `RectilinearGrid` (`loop_common.supports._3d_rectilinear_grid`). Read the golden test `test_structured_support_reference.py` for the constructor calls. + 3D grids: `G1` origin `(0,0,0)`, nsteps `(11,6,4)`, step `(1,2,3)`; `G2` origin `(100,-20,5)`, nsteps `(5,5,5)`, step `(0.5,0.5,0.5)`; `G3` origin `(-1e4,2e3,-30)`, nsteps `(21,3,8)`, step `(10,100,7.5)`. 2D: `H1` origin `(0,0)`, nsteps `(5,5)`, step `(1,1)`; `H2` origin `(10,-5)`, nsteps `(8,4)`, step `(0.5,2)`. + `nsteps` counts nodes. Extent `L = nsteps_cells * step_vector`. Length tolerance `1e-12 * max(1, L)`; counts must be exactly equal. + `test_reapplying_same_value_10_times_changes_nothing_3d`: 3 grids x 4 setters (origin, maximum, nsteps, step_vector) = 12 cells; after 10 repeats `origin`, `maximum`, `step_vector`, `nsteps` equal the start (SC-011). + `test_nsteps_setter_keeps_extent`: 3 cells; nsteps doubled then halved; `maximum` equals the start; node count equals the request. `test_step_vector_setter_encloses_extent`: 3 cells; step x 1.5; new `maximum` >= old and < old + one new step. + `test_step_vector_a_b_a_extent_drift`: 3 cells; step A, then B = 1.5 A, then A; final `maximum` minus start is >= 0 and < `A + B` on each axis; print the measured value. `test_nsteps_below_two_is_rejected_without_change`: 1 cell; `ValueError`; state equal. + `test_nsteps_not_whole_is_rejected_without_change`: 1 cell; `(5.5, 5, 5)` and NaN raise `ValueError`, state equal; `(5.0, 5.0, 5.0)` is accepted and stored as `int`. + `test_set_nelements_keeps_extent`: 3 cells; after `set_nelements(200)` `maximum` equals the start; the return value equals `n_elements`; `nsteps == max(ceil(L / h), 2)` with `h = (volume / 200) ** (1 / 3)` computed in the test. `test_set_nelements_small_request_is_clamped`: `G1.set_nelements(1)` raises nothing; every nsteps >= 2; extent kept. + `test_2d_counts_are_derived`: 2 grids x 3 checks = 6 cells; after `nsteps` assignment (a float list), `nsteps_cells`, `n_cell_x`, `n_cell_y` and `n_elements` agree with the new `nsteps`; `nsteps` is an int array; `maximum == origin + nsteps * step_vector` (legacy value). `test_rectilinear_resize_is_rejected`: 3 setters; `NotImplementedError`; nodes, origin and counts unchanged. + `StructuredGrid2DGeometry.n_elements` exists (`nsteps_cells[0] * nsteps_cells[1]`). Rule: tests must show a wrong answer (principle V). Write in ASD-STE100 Simplified Technical English. + +- [x] T070 [P] [US4] [M] Test: cell-centre sampling agrees with the grid cells for rotated boxes (U4.3; FR-015; SC-011). + - Files: `packages/loop_common/tests/test_cell_centre_sampling.py` + - Needs: T068 + - Wave: 32 + - Done when: `pytest packages/loop_common/tests/test_cell_centre_sampling.py` passes on the current tree (60 cells, "coverage, no defect", OP12). If a cell fails, record the failing ids in the report and keep the test. + - Context: New file. Boxes: `unit` and `offset_anisotropic` from `test_bounding_box_structured_grid.py` (read it for the construction). Frames: rotation 0, 30 and 90 degrees about z; local origin `(0,0,0)` and `(7,-3,11)` (6 frames; with 2 boxes = 12 cells for each test). + Set the frame with `box.set_local_transform(local_origin=..., local_rotation=...)`. Read `test_rigid_transform_contract.py` for the call and for the tolerance `Tp`. Rotation 0 with origin 0 is the control. + `test_each_cell_has_exactly_one_centre` (12 cells; the main evidence for FR-015 and SC-011): `grid = box.structured_grid(local_coordinates=True)`; `world = box.reproject(box.cell_centres(order="F"))`; `grid.cell_index_at(box.project(world))` equals `arange(n_cells)`: no `-1`, no repeat; `np.bincount` is all 1. + `test_local_coordinates_grid_associates_own_centres` (12 cells): every `box.project(world)` falls in a cell of that grid (`cell_index_at(...) != -1`). `test_centres_are_cell_midpoints_in_world` (12): `world` equals the mean of the 8 corner nodes of each cell of `box.structured_grid()` within `Tp`. + `test_centres_match_exported_cells` (12): world centres equal `grid.cell_centres` within `Tp`. `test_centres_inside_box` (12): each world centre is inside `[origin, maximum]` within `Tp`. + Only the first two tests are evidence for SC-011. Say so in their docstrings. Do not change `_bounding_box.py`. Write in ASD-STE100 Simplified Technical English. + +- [x] T071 [P] [US4] [M] Test: the pure residual and exclusion functions and the report class (U4.4, U4.5; FR-023, FR-024). + - Files: `packages/loop_interpolation/tests/test_residual_exclusion_report.py` + - Needs: T068 + - Wave: 32 + - Done when: `pytest packages/loop_interpolation/tests/test_residual_exclusion_report.py` fails before T073 with `ImportError` (the names do not exist yet), then passes after T073. + - Context: New file. No model. Import from `loop_interpolation._diagnostics`: `ConstraintDiagnosticsReport`, `ConstraintFamilyDiagnostics`, `RegionCoverageDiagnostics` (read the file for their fields and for how to build a family), and the new `direction_cosines`, `build_exclusion_records`, `summarise_residuals`, `ExclusionRecord`, `ResidualSummary`, `UNKNOWN_EXPLANATION`. + Signatures: `direction_cosines(model_vectors, data_vectors) -> ndarray (N,)`; `build_exclusion_records(excluded_rows: dict[str, dict[str, int]], families: dict[str, ConstraintFamilyDiagnostics], rows_per_point: dict[str, int] | None = None) -> tuple[ExclusionRecord, ...]`. + `test_direction_cosines_pure_function` (6 cells): equal, opposite and orthogonal vectors give 1, -1 and 0; a zero row and a NaN row give NaN; shape mismatch raises `ValueError`; both inputs unchanged; scale invariance. + `test_report_without_verification_is_unchanged` (1): `to_dict().keys()` equals the pre-S5 key set (no existing test has it: write this literal: report keys `interpolator_type`, `total_rows`, `families`, `region_coverage`, `outside_model_points`; family keys `active`, `row_count`, `dropped_rows`, `effective_weight_mean`, `effective_weight_min`, `effective_weight_max`, `source_point_count`, `outside_model_point_count`; region keys `total_support_nodes`, `active_region_nodes`, `inactive_region_nodes`, `active_fraction`; build the report with a family and a region coverage) and the old fields keep their order (`dataclasses.fields`). + `test_builder_records_keep_their_own_key` (1): `{"restricted_orthogonal": {"invalid_reference_direction": 3}}` gives a record with `family == "restricted_orthogonal"`, count 3, no `unattributed` record, no family count changed. `test_unknown_reason_keeps_its_count` (1): reason `"x"` keeps its count and gets `UNKNOWN_EXPLANATION`. + `test_unattributed_rows_are_reported` (1): `value` family, `dropped_rows` 5, outside points 2: `unattributed` count 3. `test_gradient_outside_points_give_no_false_unattributed` (1): `gradient` family, 10 source points, 2 outside, `dropped_rows` 4, `rows_per_point={"gradient": 2}`: no record; with `dropped_rows` 6: `unattributed` 2. + `test_dropped_rows_none_gives_no_unattributed` (1): `interface` family, `dropped_rows=None`, no outside points: no record. `test_surfe_style_dropped_equals_outside` (1): `rows_per_point=None`, `dropped_rows` equal to the outside count: no record. Write in ASD-STE100 Simplified Technical English. + +- [x] T072 [US4] [M] Record the S5 pre-fix evidence for the resize and cell-centre tests in `verification.md`. + - Files: `specs/005-coordinate-frame-correctness/verification.md` + - Needs: T069, T070 + - Wave: 33 + - Done when: `verification.md` has the heading `## S5 pre-fix failures` with the failing test ids, the first error of each, and the measured values; no S5 fix commit exists yet. + - Context: Run the two `loop_common` test files on the unfixed code, in one call: `pytest packages/loop_common/tests/test_resize_paths.py packages/loop_common/tests/test_cell_centre_sampling.py` (PowerShell; the worktree on `PYTHONPATH`; the venv python is `C:\Users\lachl\dev\LoopStructural\.venv\Scripts\python.exe`). + Record: the failing resize tests and the first error of each; the measured extent drift of the `nsteps` setter (example: nodes 11 to 21 with step 1 gives `maximum` 10.48 instead of 10); the extent drift of the `step_vector` setter; the stale 2D counts. Record the cell-centre tests as "coverage, no defect" if they pass, or "reproduced" with ids if not (OP12). + The 2D `maximum` is unchanged, so no 2D `maximum` test must fail. Use the layout of `## S4 pre-fix failures` in the same file. FR-031. Write in ASD-STE100 Simplified Technical English. + +## Phase 19: Stage S5 diagnostics types and resize code + +- [x] T073 [P] [US4] [M] Add the residual and exclusion types and functions to `_diagnostics.py`. + - Files: `packages/loop_interpolation/src/loop_interpolation/_diagnostics.py` + - Needs: T071 + - Wave: 33 + - Done when: `pytest packages/loop_interpolation/tests/test_residual_exclusion_report.py` passes, and `pytest packages/loop_interpolation/tests -k diagnostics` still passes. + - Context: The file has `from __future__ import annotations`; keep the existing classes and field order. Add `ResidualSummary` (frozen dataclass: `name: str`, `kind: str` ("physical" or "assembled"), `unit: str` ("field", "degree" or "assembled"), `count: int` (finite rows), `nan_count: int`, `rms`, `max_abs`, `mean_abs`: `float | None`, `None` when `count == 0`). + Add `ExclusionRecord` (frozen: `family: str`, `reason: str`, `count: int`, `explanation: str`). Add `EXCLUSION_EXPLANATIONS: dict[str, str]` with the keys `invalid_reference_value`, `invalid_reference_direction`, `evaluation_error`, `outside_model`, `unattributed` (plain one-sentence texts; `outside_model` says the count is in points, `unattributed` says it is in matrix rows) and `UNKNOWN_EXPLANATION = "No explanation is registered for this reason."`. + `summarise_residuals(name, kind, unit, residuals) -> ResidualSummary`: flatten; count non-finite values in `nan_count`; statistics from finite values only; do not change the input. + `direction_cosines(model_vectors, data_vectors) -> np.ndarray`: shape mismatch raises `ValueError`; per row dot / (norm * norm); a row with a zero or not finite norm gives NaN; inputs unchanged. + `build_exclusion_records(excluded_rows, families, rows_per_point=None) -> tuple[ExclusionRecord, ...]`, sorted by `(family, reason)`: one record per reason in `excluded_rows` (family = its key; subtract nothing); one `outside_model` record per family with `outside_model_point_count > 0`; one `unattributed` record per family with `dropped_rows` not `None` and `dropped_rows - outside_model_point_count * k > 0`, with `k = rows_per_point.get(family, 1)`; the count is the difference. Count 0 makes no record. An unknown reason gets `UNKNOWN_EXPLANATION`. + Report class: append after `outside_model_points`: `physical_residuals: dict[str, ResidualSummary] = field(default_factory=dict)`, `assembled_residuals: dict[str, ResidualSummary] = field(default_factory=dict)`, `exclusions: tuple[ExclusionRecord, ...] = ()`. `to_dict` adds `"physical_residuals"`, `"assembled_residuals"` (dicts of the fields of each summary) and `"exclusions"` (list of dicts) only when not empty. `summary()` adds one line per non-empty block after the existing lines. `total_rows` and `active_families` do not change. numpy docstrings and type hints (ruff `D`, `ANN`); py39. + +- [x] T074 [P] [US4] [S] Export four new names from `loop_interpolation`. + - Files: `packages/loop_interpolation/src/loop_interpolation/__init__.py` + - Needs: T073 + - Wave: 34 + - Done when: `python -c "from loop_interpolation import ResidualSummary, ExclusionRecord, direction_cosines, build_exclusion_records"` prints no error. + - Context: Read the file: it already imports from `._diagnostics` and has an `__all__`. Add `ResidualSummary`, `ExclusionRecord`, `direction_cosines` and `build_exclusion_records` to the same import and to `__all__`, in the style of the neighbours. Do NOT export `summarise_residuals`, `UNKNOWN_EXPLANATION` or `EXCLUSION_EXPLANATIONS`. Change no other line. + +- [x] T075 [P] [US4] [M] Keep the extent in the 3D `nsteps` and `step_vector` setters and in `set_nelements`. + - Files: `packages/loop_common/src/loop_common/geometry/_structured_grid_3d.py` + - Needs: T069, T072 + - Wave: 34 + - Done when: `pytest packages/loop_common/tests/test_resize_paths.py -k "3d or nsteps or step_vector or set_nelements or a_b_a"` passes. + - Context: Add `from typing import Optional` (no `from __future__ import annotations`; py39). `_cells_for_length(self, length, step_vector: Optional[np.ndarray] = None)`: `step_vector` defaults to the current step; the S1 rule is unchanged (`ceil(length / step - 1e-9)`; 0 cells become 2; a negative count raises `ValueError("Cannot resize the grid.")`). + `set_nelements`: `nsteps = np.maximum(np.ceil((self.maximum - self.origin) / step_vector).astype(int), 2)`. `nsteps` setter: `nsteps = np.broadcast_to(np.asarray(nsteps), self._nsteps.shape)`; if not all finite or any `nsteps != np.round(nsteps)` raise `ValueError("Cannot resize the grid. nsteps must be whole numbers.")`; `nsteps = nsteps.astype(int)`; if any `< 2` raise `ValueError("Cannot resize the grid. nsteps must be at least 2.")`; + `extent = self.nsteps_cells * self._step_vector`; `self._step_vector = extent / (nsteps - 1)`; `self._nsteps = np.array(nsteps)`; `self._on_geometry_change()`. No state changes before a raise. + `step_vector` setter: `step_vector = np.array(step_vector, dtype=float)`; `extent = self.nsteps_cells * self._step_vector`; `cells = self._cells_for_length(extent, step_vector=step_vector)`; `self._nsteps = cells + 1`; `self._step_vector = step_vector`; `self._on_geometry_change()`. + Before the edit run `git grep -n "\.nsteps = \|\.step_vector = \|set_nelements"` and read each writer; report one that depends on the old product rule. Keep the `origin` and `maximum` setters of S1. Principle: change no public name. + +- [x] T076 [P] [US4] [M] Derive the 2D counts from `nsteps` and keep the legacy `maximum`. + - Files: `packages/loop_common/src/loop_common/geometry/_structured_grid_2d.py` + - Needs: T069, T072 + - Wave: 34 + - Done when: `pytest packages/loop_common/tests/test_resize_paths.py -k 2d_counts` passes, and `pytest packages/loop_common/tests/test_2d_discrete_support.py` and `pytest tests/unit/interpolator/test_2d_discrete_support.py` still pass with no change (the four `maximum` assertions; separate calls). + - Context: First run `git grep -n "n_cell_x\|n_cell_y\|nsteps_cells\|\.maximum ="` and read each writer outside `__init__`; a writer must move to the derived form (report it). In `__init__` (lines 23 to 30) replace the assignments of `self.nsteps`, `self.maximum`, `self.nsteps_cells`, `self.n_cell_x` and `self.n_cell_y` with `self.nsteps = nsteps` (the setter coerces); keep the other lines. `nsteps` is set before `origin` is read. + Add: `nsteps` property returning `self._nsteps`; setter `self._nsteps = np.ceil(np.asarray(nsteps)).astype(int)`; read-only properties `nsteps_cells` (`self.nsteps - 1`), `n_cell_x` (`self.nsteps[0] - 1`), `n_cell_y` (`self.nsteps[1] - 1`); read-only `maximum` returning `self.origin + self.nsteps * self.step_vector`. + Keep a comment on `maximum`: legacy value (FR-028); it counts nodes as cells, so it is one step beyond the last node; `nodes` uses `nsteps_cells`; do not change it. Keep the NumPy integer type of the old attributes. Use numpy docstrings and type hints. Rule: the four existing `maximum` assertions pass with no change (decision Q1). + +- [x] T077 [P] [US4] [S] Make the resize setters of `RectilinearGrid` raise. + - Files: `packages/loop_common/src/loop_common/supports/_3d_rectilinear_grid.py` + - Needs: T069, T072 + - Wave: 34 + - Done when: `pytest packages/loop_common/tests/test_resize_paths.py -k rectilinear` passes, and `pytest packages/loop_common/tests/test_rectilinear_grid.py` still passes. + - Context: `RectilinearGrid` inherits the `origin`, `nsteps` and `step_vector` setters of the 3D geometry; a call changes the counts but not the node arrays. Add a module constant `_RESIZE_MESSAGE = "A rectilinear grid has fixed node arrays and cannot be resized. Build a new RectilinearGrid."`. + Override the three properties: each getter returns the private attribute (`self._origin`, `self._nsteps`, `self._step_vector`); each setter does `raise NotImplementedError(_RESIZE_MESSAGE)` and changes no state. Pattern: `@property def origin(self): return self._origin` then `@origin.setter def origin(self, origin): raise NotImplementedError(_RESIZE_MESSAGE)`. The constructors write the private attributes, not the setters, so construction is unchanged. `maximum` has no setter: leave it. Use numpy docstrings and type hints. + +## Phase 20: Stage S5 golden review and feature tests + +- [x] T078 [US4] [M] Survey the numeric tests, review the golden changes and regenerate the 3D golden data (C-004). + - Files: `specs/005-coordinate-frame-correctness/verification.md`, `packages/loop_common/tests/data/structured_support_reference.npz`, `packages/loop_common/tests/data/structured_grid_2d_reference.npz` + - Needs: T075, T076, T077 + - Wave: 35 + - Done when: `pytest packages/loop_common/tests/test_structured_support_reference.py packages/loop_common/tests/test_structured_grid_2d_reference.py` passes, and `verification.md` has the heading `## S5 golden review (C-004)` with one entry for each changed key or test (old value, new value, cause). + - Context: Order (C-004, SC-015: 0 tolerance changes). (1) Run the 3D golden test and list each failing key. (2) Check each key against the rule: each setter keeps the extent or encloses it by less than one cell. (3) Write each key with its old and new value in the section. (4) Only then run `python packages/loop_common/tests/test_structured_support_reference.py --update`. (5) Run the test again. + Run the 2D golden test. It must not change (the 2D `maximum` stays). Run `--update` for it only if a key differs, and record the key. Numeric survey, in separate calls: `pytest tests/unit/modelling`; `pytest packages/loop_interpolation/tests`; `pytest packages/loop_common/tests` (PowerShell; worktree on `PYTHONPATH`; `packages/loopstructural_visualisation/src` also on `PYTHONPATH`; `pyevtk` is only in a scratch folder: add it to `PYTHONPATH` or report the skip). + For each test that fails or changes a number, record the test, the old value, the new value and the cause. Do not edit a test file. `packages/loop_interpolation/tests/golden` is the parity oracle (OP17): record each change there and do NOT regenerate it. If the cause is the intended extent change, name the file for the maintainer; if the cause is not clear, report it. Write in ASD-STE100 Simplified Technical English. + +- [x] T079 [P] [US4] [M] Record the S5 post-fix evidence for the resize and cell-centre tests. + - Files: `specs/005-coordinate-frame-correctness/verification.md` + - Needs: T078 + - Wave: 36 + - Done when: `verification.md` has the heading `## S5 resize and cell-centre results` and `pytest packages/loop_common/tests/test_resize_paths.py packages/loop_common/tests/test_cell_centre_sampling.py` passes. + - Context: Run the two files in one call. Record the commands, the counts (resize: 12 + 3 + 3 + 3 + 1 + 1 + 3 + 1 + 6 + 3 = 36 cells; cell centres: 5 tests x 12 = 60 cells), the times, and the measured extent drift of the A, B, A test for each grid (the bound is `A + B`). Compare with `## S5 pre-fix failures`: for each id that failed, write that it now passes. + Write that `test_each_cell_has_exactly_one_centre` and `test_local_coordinates_grid_associates_own_centres` are the evidence for FR-015 and SC-011, and the other three are coverage. Record OP12: "not reproduced" if all 60 cells passed before the fix, and no change to `_bounding_box.py`. Record the `ruff check LoopStructural packages` result for the three `loop_common` source files. Write in ASD-STE100 Simplified Technical English. + +- [x] T080 [US4] [M] Test: misfit methods keep state, compare in world coordinates, and keep their shapes (U4.4; FR-023; SC-009). + - Files: `tests/unit/modelling/test_diagnostics_frames.py` + - Needs: T073, T078 + - Wave: 36 + - Done when: `pytest tests/unit/modelling/test_diagnostics_frames.py -k "leaves_state_unchanged or exact_field or value_ or empty or old_misfit_shapes"` fails before T083 (record the failing ids and the first error of each) and passes after T084. + - Context: New file. Use `world_foliation_model`, `ANCHORS`, `ROTATIONS_DEG` from `tests.fixtures.coordinate_frames`; solve with the pattern of `tests/unit/modelling/test_reference_constraint_frames.py` (read both). Engines: FDI and PLI. Cell grid: 4 anchors x 3 rotations x 2 engines = 24 (value tests: 12 cells). + `test_misfit_leaves_state_unchanged` (24): after a solve, take deep copies of every `interpolator.data[...]` array, `c`, `interpolator.up_to_date`, `builder._up_to_date`, the weights of `constraints` and `id(latest_diagnostics_report)`; call the two OLD methods `evaluate_gradient_misfit`, `evaluate_value_misfit`; every snapshot is equal (`np.array_equal`). It uses no new name. `test_new_misfit_methods_leave_state_unchanged` (24): the same, with `evaluate_gradient_cosines`, `evaluate_value_residuals` and `verification_report()`. + `test_gradient_cosine_matches_analytic_direction_exact_field` (24): `monkeypatch` replaces `feature.evaluate_gradient` with a stub that returns the analytic world direction of the R1 foliation at the given points; `abs(1 - cos) <= 1e-10` for the old method (one row) and for `evaluate_gradient_cosines`. `..._solved_model` (24): the real solved model; `tol_solved` is measured: record the maximum and set the bound with a margin of 2 (it is not `1e-10`). + `test_value_residuals_are_absolute_difference` (12): equals `abs(value - evaluate_value(world_points))` at the 27 value rows; no NaN. `test_value_misfit_is_normalised_by_field_range` (12): equals `evaluate_value_residuals() / (feature.max() - feature.min())` within `1e-12`. `test_value_misfit_with_zero_range_does_not_divide` (2): `monkeypatch` makes `feature.max` and `feature.min` return the same number (cell 1) or NaN (cell 2); `warnings.simplefilter("error")`; the result equals the unnormalised absolute difference. + `test_misfit_with_no_rows_is_empty` (2): no gradient rows gives shape `(0,)`; no value rows gives shape `(0,)`. `test_old_misfit_shapes_are_kept` (4): one family `(1, N)`; two families of equal length `(2, N)`; two of different length `(1, N1 + N2)` (OP16, default); values equal `evaluate_gradient_cosines()`. + +- [x] T081 [US4] [M] Test: a stale model is not rebuilt, and the report keeps residuals and exclusions apart (U4.4, U4.5; FR-023, FR-024; SC-009). + - Files: `tests/unit/modelling/test_diagnostics_frames.py` + - Needs: T080 + - Wave: 37 + - Done when: `pytest tests/unit/modelling/test_diagnostics_frames.py -k "stale or report_blocks or unsolved or exclusions or outside_points"` fails before T083 and T084 (record the failing ids; `test_stale_misfit_is_not_rebuilt` must fail on the current tree because the call rebuilds the model) and passes after T084. + - Context: Append to the file of T080 (same fixtures; engines FDI and PLI). Stale recipe (settled, plan "Answers"): after the solve, for a stale interpolator assign `interpolator.up_to_date = False`; for a stale builder assign `builder._up_to_date = False` directly (not `builder.set_not_up_to_date`: it also stales the interpolator); assert the flag before the call. PLI note (settled): `P1Interpolator.add_gradient_constraints` is `pass`, so PLI assembles no `gradient` rows and the builder reports `dropped_rows = 2 * source_points` for it; `unattributed` for PLI gradient data is true. Run the no-false-unattributed check with FDI only; for PLI record in the test output the measured `unattributed` count and the `outside_model` count (print them) and assert `outside_model == 3`. `test_stale_misfit_is_not_rebuilt` (2 engines x 2 kinds, interpolator stale or builder stale = 4): uses the OLD methods only. Snapshot `interpolator.up_to_date`, `builder._up_to_date`, `c` and the data arrays; after the call every snapshot is equal (on the current tree `up_to_date` changes from `False` to `True`: the test fails; FR-031). + `test_stale_misfit_warns_once_and_returns_nan` (4): the old and the new methods each write exactly one `UserWarning` that names the stale object; the new methods return `np.full(n, np.nan)`; the old methods return the old shape filled with NaN. `test_stale_verification_report_does_not_rebuild` (4): snapshot equal; both residual blocks empty; exclusions present; one warning. + `test_report_blocks_are_separate` (2): the two dicts are different objects; physical entries have `kind == "physical"`, assembled entries `kind == "assembled"`; the physical `value` entry has unit `"field"` and its `rms` equals the RMS of `evaluate_value_residuals()` within `1e-12`; the assembled entries equal `summarise_residuals(name, "assembled", "assembled", r)` for `calculate_residual_for_constraints()`; `to_dict` has two separate blocks. Import `summarise_residuals` and `UNKNOWN_EXPLANATION` from `loop_interpolation._diagnostics`. + `test_unsolved_model_gives_no_residual_blocks` (2): `interpolator.c = None`; both blocks `{}`; exclusions present; no exception. `test_exclusions_are_counted_and_explained` (2): the restricted-support reference case of `test_reference_constraint_frames.py`: builder records under `restricted_equality` and `restricted_orthogonal` with the analytic count; the explanation is not empty and not `UNKNOWN_EXPLANATION`. + `test_gradient_outside_points_give_no_false_unattributed_model` (2): gradient data with 3 points outside the box: an `outside_model` record of 3 and no `unattributed` record for `gradient`. + +- [x] T082 [US4] [M] Record the S5 pre-fix evidence for the diagnostics tests in `verification.md`. + - Files: `specs/005-coordinate-frame-correctness/verification.md` + - Needs: T081 + - Wave: 38 + - Done when: `verification.md` has the heading `## S5 diagnostics pre-fix failures` with the failing ids, the first error of each, and the measured values; no misfit fix commit exists yet. + - Context: Run `pytest tests/unit/modelling/test_diagnostics_frames.py` on the code before T083 (root suite only; one call). Record: the state test with the old methods (a different array in `data["gradient"]` after the call); the exact-field cosine test (one column, local vectors); the stale-state test (the call rebuilds the model: `up_to_date` goes from `False` to `True`); the tests of the new names (they fail because the names do not exist). + Record `test_value_residuals_are_absolute_difference` as "not reproduced" if it passes (the value path does not change the stored positions). Add one behavioural measure with the production path: for `ANCHORS[1]` and rotation 30, the largest absolute change of `interpolator.data["gradient"]` after one call of `evaluate_gradient_misfit`. Use the layout of `## S4 pre-fix failures`. FR-031. Write in ASD-STE100 Simplified Technical English. + +## Phase 21: Stage S5 misfit code + +- [x] T083 [US4] [M] Rewrite the misfit methods on shared non-mutating code and add the new methods. + - Files: `LoopStructural/modelling/features/_geological_feature.py` + - Needs: T073, T082 + - Wave: 39 + - Done when: `pytest tests/unit/modelling/test_diagnostics_frames.py -k "not report_blocks and not unsolved and not exclusions and not outside_points and not verification_report and not new_misfit_methods"` passes. + - Context: The module has `from __future__ import annotations`. Add (import `direction_cosines` from `loop_interpolation._diagnostics`; `warnings`): `_stale_reason(self) -> Optional[str]`: text for the first stale object among `self.builder._up_to_date` is False, each builder in `self.builder.faults`, `self.interpolator.up_to_date` is False; `None` otherwise (builder `None`: only the interpolator). The text names the object and says to build the model first. + `_gradient_cosines_fresh(self) -> list`: for each of `get_gradient_constraints()` and `get_norm_constraints()` with rows: `rows = np.array(array, dtype=float, copy=True)`; `positions = rows[:, :3]`, reproject with `self.interpolator.bounding_box.reproject` when `self.interpolator.bounding_box` is not `None` (settled: `GeologicalFeature` has no `bounding_box`; do not use `self.model.bounding_box`); `data = rows[:, 3:6]`, convert with `self.interpolator.local_to_world_vectors(data)` when that bounding box is not `None`; `model = self.evaluate_gradient(positions)`; append `direction_cosines(model, data)`. + `_value_residuals_fresh(self) -> np.ndarray`: copy `get_value_constraints()`, reproject positions, return `np.abs(rows[:, 3] - self.evaluate_value(positions))`. + `evaluate_gradient_cosines() -> ndarray (N,)`: on a stale model `warnings.warn(text, UserWarning, stacklevel=2)` and return `np.full(n, np.nan)` (`n` = row count of the stored arrays); else `np.concatenate` of the list (`np.zeros(0)` if empty). `evaluate_value_residuals() -> ndarray (N,)`: same rule with `_value_residuals_fresh`. + `evaluate_gradient_misfit()` (old name, new body): same code; one row per non-empty family: `(1, N)` for one family, `(2, N)` for two of equal length, `(1, N1 + N2)` for two of different length, `(0,)` with no rows; stale gives that shape filled with NaN. `evaluate_value_misfit()`: stale gives NaN `(N,)`; else the residuals divided by `self.max() - self.min()` only when the range is finite and > 0, else unnormalised. + Remove every `self.builder.up_to_date()` call and the in-place `grad /= norm` from these methods. Never write to a stored array, weight, coefficient or flag. numpy docstrings (state that the gradient value is a cosine in world coordinates; see COMPAT.md) and type hints. Read the whole old methods first (about lines 240 to 300). + +- [x] T084 [US4] [M] Add `GeologicalFeature.verification_report` with physical and assembled residuals and exclusions. + - Files: `LoopStructural/modelling/features/_geological_feature.py` + - Needs: T083 + - Wave: 40 + - Done when: `pytest tests/unit/modelling/test_diagnostics_frames.py` passes (all cells), and `pytest tests/unit/modelling/test_scalar_field_grid.py` still passes. + - Context: Add `_physical_residual_summaries(self) -> dict[str, ResidualSummary]`: `value` from `_value_residuals_fresh()` (kind `"physical"`, unit `"field"`); `gradient` and `normal` from `_gradient_cosines_fresh()`: `np.degrees(np.arccos(np.clip(cosine, -1, 1)))`, unit `"degree"`; use `summarise_residuals(name, "physical", unit, values)`. Add the public `verification_report(self) -> ConstraintDiagnosticsReport`. + Body: `base = self.interpolator._build_constraint_diagnostics_report()` (a new object; do NOT set `latest_diagnostics_report`). `physical = {}`, `assembled = {}`. If `_stale_reason()` is not `None`: `warnings.warn(text, UserWarning, stacklevel=2)` and skip both blocks. Else `physical = self._physical_residual_summaries()`; and only if `interpolator.c` is not `None` and the interpolator has `calculate_residual_for_constraints`: for `(name, residual)` in `calculate_residual_for_constraints().items()` set `assembled[name] = summarise_residuals(name, "assembled", "assembled", residual)`. + Exclusions are always built (they read state only): `excluded = getattr(self.builder, "excluded_rows", {})`; `rows_per_point = {"gradient": 2}` if the interpolator is a `DiscreteInterpolator`, else `None`. Return `dataclasses.replace(base, physical_residuals=physical, assembled_residuals=assembled, exclusions=build_exclusion_records(excluded, base.families, rows_per_point))`. + An unsolved model (`c` is `None`) leaves BOTH residual blocks empty with no warning (plan decision 16; settled after the wave 37 review: guard the physical block on `interpolator.c is not None` too, because `evaluate_value` with `c = None` fails; the test `test_unsolved_model_gives_no_residual_blocks` requires both `{}`). The stale warning text must contain the word `interpolator` or `builder` (the stale object). The `gradient` rows_per_point is 2 for any `DiscreteInterpolator` (FDI gives 2 rows for each point; PLI gives 0 rows, so PLI gradient data gives an `unattributed` record, measured 8). Do not change `calculate_residual_for_constraints`. numpy docstring and type hints. Import `ConstraintDiagnosticsReport`, `ResidualSummary`, `summarise_residuals`, `build_exclusion_records` from `loop_interpolation._diagnostics`. + +## Phase 22: Stage S5 finish + +- [x] T085 [P] [Finish] [M] Record the S5 notices in `COMPAT.md`. + - Files: `COMPAT.md` + - Needs: T084 + - Wave: 41 + - Done when: `grep -c "evaluate_gradient_cosines" COMPAT.md` prints at least 1 and `grep -c "RectilinearGrid" COMPAT.md` prints at least 1. + - Context: Read `COMPAT.md` first, the S1 to S4 notices, and copy their layout. Add migration notices. (1) `GeologicalFeature.evaluate_gradient_misfit()` changed the stored constraint arrays in place, compared a local vector with a world vector, and read one column. It keeps its name and its public shape and changes nothing. It returns the cosine between the world data direction and the model gradient for each row. The old value was a dot product with an unnormalised model gradient, so the new value differs. Use `evaluate_gradient_cosines()` for a 1-D array. + Two families of different length now give one row `(1, N1 + N2)`; before, this raised `ValueError` (OP16). (2) Both misfit methods no longer rebuild a stale model: they write one `UserWarning` and return NaN; build the model first. (3) `evaluate_value_misfit()` keeps its normalisation by the field range and no longer divides by zero (zero or not finite range gives the absolute difference). `evaluate_value_residuals()` gives the unnormalised 1-D result. + (4) Assigning `nsteps` or `step_vector` on a 3D structured grid, and `set_nelements`, keep the extent (a new step, not a product rule); `nsteps` must be whole and at least 2. `set_nelements` keeps `maximum`. Grid steps and solved values can change slightly. (5) A `RectilinearGrid` raises `NotImplementedError` on assignment of `origin`, `nsteps` or `step_vector`. (6) The 2D grid derives `nsteps_cells`, `n_cell_x`, `n_cell_y` from `nsteps`; `maximum` keeps its legacy value. (7) `ConstraintDiagnosticsReport` has three new optional fields (`physical_residuals`, `assembled_residuals`, `exclusions`); with an identity frame and no resize the value misfit is the same as before. + Do not change other text. Write in ASD-STE100 Simplified Technical English. + +- [x] T086 [P] [Finish] [S] Record the S5 misfit decision in `ROADMAP.md`. + - Files: `ROADMAP.md` + - Needs: T084 + - Wave: 41 + - Done when: `grep -c "Misfit diagnostics of GeologicalFeature (spec 005, S5)" ROADMAP.md` prints 1. + - Context: Add a new `###` subsection as the last one of `## Decisions`, directly before `## Stage sequence`. Heading: `### Misfit diagnostics of GeologicalFeature (spec 005, S5)`. Body (ASD-STE100): decision date 2026-11-05 (spec 005, S5, Q2); the misfit methods `evaluate_gradient_misfit` and `evaluate_value_misfit` keep their names and public shape; the gradient value is now a cosine in world coordinates; a stale model gives a `UserWarning` and NaN (no rebuild); + the new names `evaluate_gradient_cosines` and `evaluate_value_residuals` give the 1-D result; the notice is in `COMPAT.md`. Keep the style of the neighbouring subsections. Do not change other text. + +- [x] T087 [Finish] [M] Record the S5 rollback check in `verification.md`. + - Files: `specs/005-coordinate-frame-correctness/verification.md` + - Needs: T074, T084 + - Wave: 41 + - Done when: `verification.md` has the heading `## S5 rollback checks` with the results of the 6 checks below, and `git worktree list` shows no throw-away worktree. + - Context: Use a throw-away worktree at the S4 end commit (the commit before the first S5 code commit): `git worktree add --detach ../LoopStructural-005-s5-rollback `, with its own `PYTHONPATH` (check `LoopStructural.__file__` and `loop_interpolation.__file__`). Checks: (a) `dataclasses.fields` of `ConstraintFamilyDiagnostics`, `RegionCoverageDiagnostics` and the first five fields of `ConstraintDiagnosticsReport` are equal on both trees; + (b) `to_dict()` of a report built with no verification is equal on both trees; (c) `from loop_interpolation import *` on the S4 tree gives a subset of the S5 names; (d) the signatures of `GeologicalFeature.evaluate_gradient_misfit` and `evaluate_value_misfit` are equal, and the result shape for one family is `(1, N)` on both; + (e) with an unchanged grid, `origin`, `maximum`, `nsteps` and `step_vector` of a structured support are equal on both trees; (f) the four 2D `maximum` assertions pass on both trees. Then `git worktree remove --force ` and `git worktree prune`. Write in ASD-STE100 Simplified Technical English. + +- [x] T088 [Finish] [M] Run the S5 gates and complete `verification.md` and the `ROADMAP.md` status log. + - Files: `specs/005-coordinate-frame-correctness/verification.md`, `ROADMAP.md` + - Needs: T079, T085, T086, T087 + - Wave: 42 + - Done when: `ruff check LoopStructural packages && pytest` and `pytest packages/loop_common/tests packages/loop_interpolation/tests` pass, `verification.md` has the S5 sections named below, and `ROADMAP.md` has one new status-log line for S5. + - Context: Run the root `pytest` and the package tests in separate calls (the conftest plugins clash) with the worktree `PYTHONPATH`, plus `packages/loopstructural_visualisation/src`; `pyevtk` is only in a scratch folder (add it or report the skip). Add to `verification.md`: `## S5 post-fix results` (commands, counts, times, the cells of `test_diagnostics_frames.py` and `test_residual_exclusion_report.py`, the measured `tol_solved`), + `## S5 backend matrix` (FDI and PLI misfit, residual and exclusion report: mandatory; resize of structured, 2D and rectilinear grids: mandatory; P2, native, Surfe, residual reports for fold and P2 cells: not in S5; S2 and S6 cover them), `## S5 test run time (OP2)`, and a line for each open point: OP12 (cell-centre result), OP16 (default shape kept), OP17 (parity files changed or not, from T078). + Add to the status log of `ROADMAP.md` a line dated today: "Spec 005 stage S5 implemented: resize paths keep the extent, cell-centre guard tests, non-mutating world-correct misfit diagnostics, residual and exclusion report." Report any failing gate and any file under `packages/loop_interpolation/tests/golden` that changed. Write in ASD-STE100 Simplified Technical English. + ## Wave table | Wave | Tasks | Tier [S] | Tier [M] | @@ -1201,10 +1416,21 @@ frame. A legacy file needs an explicit choice. A frame change never gives a stal | 29 | T063, T064, T065, T066 | 2 | 2 | | 30 | T067 | 0 | 1 | | 31 | T068 | 0 | 1 | +| 32 | T069, T070, T071 | 0 | 3 | +| 33 | T072, T073 | 0 | 2 | +| 34 | T074, T075, T076, T077 | 2 | 2 | +| 35 | T078 | 0 | 1 | +| 36 | T079, T080 | 0 | 2 | +| 37 | T081 | 0 | 1 | +| 38 | T082 | 0 | 1 | +| 39 | T083 | 0 | 1 | +| 40 | T084 | 0 | 1 | +| 41 | T085, T086, T087 | 1 | 2 | +| 42 | T088 | 0 | 1 | ## Critical path -T003 -> T004 -> T005 -> T010 -> T013 -> T015 -> T017 -> T018 -> T020 -> T023 -> T025 -> T030 -> T031 -> T033 -> T035 -> T036 -> T039 -> T040 -> T041 -> T043 -> T047 -> T050 -> T051 -> T052 -> T053 -> T054 -> T055 -> T057 -> T063 -> T067 -> T068 +T003 -> T004 -> T005 -> T010 -> T013 -> T015 -> T017 -> T018 -> T020 -> T023 -> T025 -> T030 -> T031 -> T033 -> T035 -> T036 -> T039 -> T040 -> T041 -> T043 -> T047 -> T050 -> T051 -> T052 -> T053 -> T054 -> T055 -> T057 -> T063 -> T067 -> T068 -> T069 -> T072 -> T075 -> T078 -> T080 -> T081 -> T082 -> T083 -> T084 -> T085 -> T088 ## Dependencies @@ -1239,4 +1465,8 @@ T003 -> T004 -> T005 -> T010 -> T013 -> T015 -> T017 -> T018 -> T020 -> T023 -> (`rescale` and `scale` warn), T060 and T061 (exporters), T062 (example). - Finish: T063 (golden review) after the code; T064 to T066 (`COMPAT.md`, `ROADMAP.md`, `API.md`); T067 (rollback check); T068 (gates) last. -- Stages S5 and S6 have no tasks yet (OP1 in plan.md). +- S5 starts after T068 (S4 is its base). Tests, wave 32: T069 (resize), T070 (cell centres), T071 (residual and exclusion functions). T072 (pre-fix evidence) after T069 and T070. +- S5 code: T073 (`_diagnostics.py`) -> T074 (export); T075 (3D setters), T076 (2D counts) and T077 (rectilinear guard) after T072. T078 (numeric survey, golden review, regenerate golden data) after the three; T079 (post-fix evidence) after T078. +- S5 feature tests, one file in a chain: T080 -> T081 (after T073 and T078); T082 (pre-fix evidence). Misfit code: T083 -> T084 (one file; after T073 and T082). +- Finish: T085 (`COMPAT.md`), T086 (`ROADMAP.md`) and T087 (rollback check) after T084; T088 (gates) last. +- No open point blocks a task. OP12, OP16 and OP17 have defaults (OP17 can block the S5 merge if a parity file changes: T078 reports it). Stage S6 has no tasks yet (OP1 in plan.md). diff --git a/specs/005-coordinate-frame-correctness/verification.md b/specs/005-coordinate-frame-correctness/verification.md index 481d3bfc..441bc11f 100644 --- a/specs/005-coordinate-frame-correctness/verification.md +++ b/specs/005-coordinate-frame-correctness/verification.md @@ -697,3 +697,358 @@ of `rescale` or `scale` appears, so no S4 code path calls them. - `tests/unit/io/test_exporters.py::test_write_vol_gocad_is_broken_for_real_bounding_box` failed only because two classes are named `LoopValueError`. `BoundingBox.__getitem__` raises the class of `loop_common.geometry._bounding_box` (a `ValueError`), and the test expected the class of `LoopStructural.utils.exceptions`. The test now accepts `(LoopValueError, ValueError)`. The file passes: 13 passed (run with `pyevtk` on `PYTHONPATH`). The full root suite was not run again after this change; the failure was the only one in the S4 root run. - The task text and the S4 notice in `COMPAT.md` named `Loop3DView.add_data`. The method is `Loop2DView.add_data`. `plan.md`, `tasks.md` and `COMPAT.md` now say so. + +## S5 pre-fix failures + +Method. Both `loop_common` test files ran on the unfixed code, in one call: `python -m pytest packages/loop_common/tests/test_resize_paths.py packages/loop_common/tests/test_cell_centre_sampling.py -p no:cacheprovider -q -rf` (PowerShell, the worktree `LoopStructural-005-s5` on `PYTHONPATH`). No S5 library change exists at this commit. The result is 24 failed and 76 passed. All 24 failures are in the resize file. All 60 cells of the cell-centre file pass. + +Resize file (`test_resize_paths.py`): 24 failed, 76 passed. + +| Test function | Failed / total | Representative id and first error | +|---|---|---| +| `test_nsteps_setter_keeps_extent` | 3 / 3 | `[G1]`: the `maximum` is [10.5, 11, 10.5] and not [10, 10, 9] after `nsteps` is set to [22, 12, 8]. `[G2]` gives 102.25, -17.75, 7.25 (start 102, -18, 7). `[G3]` gives -9795, 2250, 26.25 (start -9800, 2200, 22.5) | +| `test_step_vector_setter_encloses_extent` | 2 / 3 | `[G1]`: the `maximum` is [10.5, 9, 9] and is below the start [10, 10, 9], so the extent is lost. `[G3]` fails in the same way (the `maximum` is below the start). `[G2]` passes | +| `test_nsteps_below_two_is_rejected_without_change` | 1 / 1 | `DID NOT RAISE ValueError`. The setter accepts nsteps [1, 6, 4] | +| `test_nsteps_not_whole_is_rejected_without_change` | 1 / 1 | `DID NOT RAISE ValueError`. The setter accepts [5.5, 5, 5] and [nan, 5, 5] | +| `test_set_nelements_keeps_extent` | 3 / 3 | `[G1]`: the `maximum` is not the start `maximum`. `set_nelements(200)` gives [9.43, 10.29, 10] and not [10, 10, 9]. `[G2]` gives 102.08, -17.92, 7.08. `[G3]` gives -9811, 2270, 10 | +| `test_set_nelements_small_request_is_clamped` | 1 / 1 | `assert np.all(nsteps >= 2)` fails: `nsteps` is [2, 2, 1], and the `maximum` is [5.5, 6, 0] | +| `test_2d_counts_are_derived` | 10 / 10 | All 5 checks fail for `[H1]` and for `[H2]`: `cells`, `n_cell`, `n_elements`, `dtype`, `maximum`. See the measured values below | +| `test_rectilinear_resize_is_rejected` | 3 / 3 | `[nsteps]`: `DID NOT RAISE NotImplementedError`. `[step_vector]` and `[origin]` fail in the same way | +| `test_reapplying_same_value_10_times_changes_nothing_3d` | 0 / 12 | Passes now (coverage, no defect) | +| `test_step_vector_a_b_a_extent_drift` | 0 / 3 | Passes now (coverage, no defect). See the measured drift below | + +Measured values (a short script, not saved in the repository). + +| Measure | Value | +|---|---| +| `nsteps` setter, nodes 11 to 21, step 1, origin 0 | `maximum` is 10.476 (10.48) in each axis, and the expected value is 10 | +| `step_vector` setter, step x 1.5 on G1 | `maximum` is [10.5, 9, 9], the start is [10, 10, 9]. The `nsteps` become [8, 4, 3] | +| `step_vector` setter, step x 1.5 on G3 | `maximum` is [-9805, 2150, 26.25], the start is [-9800, 2200, 22.5]. The `nsteps` become [14, 2, 6] | +| `step_vector` setter, step x 1.5 on G2 | `maximum` is [102.25, -17.75, 7.25], the start is [102, -18, 7]. The extent is enclosed, so the test passes | +| A-B-A `step_vector` drift (A, then 1.5 A, then A) | G1 [1, 0, 3]. G2 [0.5, 0.5, 0.5]. G3 [0, 0, 7.5]. All are inside the limit (A + B), so the test passes | +| 2D `nsteps` set to [nx + 3.0, ny + 2.0], H1 (5, 5) | `nsteps_cells` stays [4, 4]. `n_cell_x` and `n_cell_y` stay 4. `n_elements` stays 16. The dtype is float64. The `maximum` stays [5, 5] and the expected value is [8, 7] | +| Same set on H2 (8, 4) | `nsteps_cells` stays [7, 3]. `n_elements` stays 21. The dtype is float64. The `maximum` stays [14, 3] and the expected value is [15.5, 7] | + +The 2D `maximum` check fails only because the counts are stale after the set. The legacy formula (origin + nsteps x step) is the expected value in the test, so the fix does not change the formula. + +Cell-centre file (`test_cell_centre_sampling.py`): 60 passed, 0 failed. Result: coverage, no defect. OP12 is not reproduced. + +Note on the one-to-one test. The first version of the test used `structured_grid(local_coordinates=True)` and failed in 8 of 12 cells (rotation 30 and 90). That grid is the axis-aligned enclosure of the projected corners, and it has a different step. The premise in the task text was wrong. The test now uses the midpoints of the exported world grid and the own cell grid of the box. `_bounding_box.py` was not changed. Also, `set_local_transform` takes the argument `rotation_matrix=`. + +## S5 golden review (C-004) + +Order (C-004): (1) run the 3D golden test and list the failing keys; (2) check each key against the rule; (3) write the keys here; (4) run `--update`; (5) run the test again. SC-015: no tolerance was changed (`rtol=1e-9`, `atol=1e-12` stay). + +Step 1. `test_structured_support_reference.py` failed in 3 cells (StructuredGrid, TetMesh, P2TetMesh). A key-by-key diff (script in the scratchpad, not in the repository) found 78 changed keys of 611. The key sets are equal. Each of the 3 classes has the same 26 changed keys: 3 grids (`unit`, `aniso`, `small`) x (`step_vector`, `to_dict.step_vector`, and `after_set_nelements`, `after_nsteps`, `after_step_vector`, `after_origin`, `after_maximum` `.step_vector`) plus `after_maximum.maximum` and `after_maximum.nsteps` where they differ (`unit` has no `nsteps` change, `aniso` and `small` do). The 2D golden (`test_structured_grid_2d_reference.py`) passes with no change (4 passed): the 2D `maximum` is the same, so no 2D key changed and `--update` was not run for it. + +Step 2. Rule: each setter keeps the extent or encloses it by less than one cell. Extent = `(nsteps - 1) x step_vector`. Checked on the `StructuredGrid` keys (the 3 classes give the same numbers): + +| Grid | Stage | Old extent | New extent | Result | +|---|---|---|---|---| +| unit | `set_nelements(200)` | 3.571 | 4.000 (start 4) | kept | +| unit | `nsteps` + 1 | 4.286 | 4.000 | kept | +| unit | `step_vector` x 2 | 4.286 | 4.000 | kept | +| unit | origin - 1 | max 4.714 | max 4.333 (1 cell = 1.333; enclosed by 0.333) | enclosed by less than one cell | +| unit | `maximum` + 3 (request 7.333) | 9.000 | 8.333 (over by 1.0, 1 cell = 1.333) | enclosed by less than one cell | +| aniso | `set_nelements(200)` | [12.1, 5.18, 6.29] | [13, 5, 6.5] (start [13, 5, 6.5]) | kept | +| aniso | `nsteps` + 1, `step_vector` x 2 | [12.8, 6.09, 7.14] | [13, 5, 6.5] | kept | +| aniso | origin - 1 | max [13.2, 6.91, 7.86] | max [13.5, 5.6, 7.0] (cells [1.5, 1.6, 1.5]) | enclosed by less than one cell | +| aniso | `maximum` + 3 (request [16.5, 8.6, 10]) | [17.4, 10.55, 11.29] | [16.5, 8.8, 10] | enclosed by less than one cell | +| small | `set_nelements(200)`, `nsteps` + 1, `step_vector` x 2 | 2.143, 2.571, 2.571 | 2.000 | kept | +| small | origin - 1 | max 3.286 | max 2.333 (1 cell = 0.667) | enclosed by less than one cell | +| small | `maximum` + 3 (request 5.333) | 6.714 | 5.667 (over by 0.333, 1 cell = 0.667) | enclosed by less than one cell | + +All new values obey the rule. The old values did not keep the extent (they drifted), and that is the S5 defect. + +Causes. + +- C1. `unit/step_vector` and `to_dict.step_vector` (and the same keys for `aniso`, `small`). Old value: the step after a later setter (for example 0.7143, not the constructor value 1.0). New value: 1.0, which is the constructor step. The old `nsteps` setter did `self._step_vector /= change_factor`, an in-place change. The golden test stores the array that the property returns, so a later setter changed the stored "start" value. The new setters give a new array. The new value is the correct start value. +- C2. `after_set_nelements`, `after_nsteps`, `after_step_vector`, `after_origin`, `after_maximum` `.step_vector`: the setters keep the extent, so the new step is the extent divided by the new cell count. +- C3. `after_maximum.maximum` and `after_maximum.nsteps`: the grid before this stage has a different step and `nsteps` (C2), so the encloses-the-request result differs. The result is enclosed by less than one cell. + +Step 3. The 78 changed keys (old value, new value, cause). + +| Key | Old value | New value | Cause | +|---|---|---|---| +| `P2TetMesh/aniso/after_maximum.maximum` | [17.4, 10.5455, 11.2857] | [16.5, 8.8, 10] | C3 (maximum request is enclosed by less than one cell) | +| `P2TetMesh/aniso/after_maximum.nsteps` | [7, 9, 7] | [6, 9, 7] | C3 (maximum request is enclosed by less than one cell) | +| `P2TetMesh/aniso/after_maximum.step_vector` | [1.4, 1.81818, 1.71429] | [1.5, 1.6, 1.5] | C2 (resize keeps the extent: new step = extent / cells) | +| `P2TetMesh/aniso/after_nsteps.step_vector` | [0.7, 0.909091, 0.857143] | [0.75, 0.8, 0.75] | C2 (resize keeps the extent: new step = extent / cells) | +| `P2TetMesh/aniso/after_origin.step_vector` | [1.4, 1.81818, 1.71429] | [1.5, 1.6, 1.5] | C2 (resize keeps the extent: new step = extent / cells) | +| `P2TetMesh/aniso/after_set_nelements.step_vector` | [0.7, 0.909091, 0.857143] | [1, 0.888889, 0.9] | C2 (resize keeps the extent: new step = extent / cells) | +| `P2TetMesh/aniso/after_step_vector.step_vector` | [1.4, 1.81818, 1.71429] | [1.5, 1.6, 1.5] | C2 (resize keeps the extent: new step = extent / cells) | +| `P2TetMesh/aniso/step_vector` | [0.7, 0.909091, 0.857143] | [0.5, 2, 1.5] | C1 (old array was changed in place by a later setter) | +| `P2TetMesh/aniso/to_dict.step_vector` | [0.7, 0.909091, 0.857143] | [0.5, 2, 1.5] | C1 (old array was changed in place by a later setter) | +| `P2TetMesh/small/after_maximum.maximum` | [6.71429, 6.71429, 6.71429] | [5.66667, 5.66667, 5.66667] | C3 (maximum request is enclosed by less than one cell) | +| `P2TetMesh/small/after_maximum.nsteps` | [10, 10, 10] | [11, 11, 11] | C3 (maximum request is enclosed by less than one cell) | +| `P2TetMesh/small/after_maximum.step_vector` | [0.857143, 0.857143, 0.857143] | [0.666667, 0.666667, 0.666667] | C2 (resize keeps the extent: new step = extent / cells) | +| `P2TetMesh/small/after_nsteps.step_vector` | [0.428571, 0.428571, 0.428571] | [0.333333, 0.333333, 0.333333] | C2 (resize keeps the extent: new step = extent / cells) | +| `P2TetMesh/small/after_origin.step_vector` | [0.857143, 0.857143, 0.857143] | [0.666667, 0.666667, 0.666667] | C2 (resize keeps the extent: new step = extent / cells) | +| `P2TetMesh/small/after_set_nelements.step_vector` | [0.428571, 0.428571, 0.428571] | [0.4, 0.4, 0.4] | C2 (resize keeps the extent: new step = extent / cells) | +| `P2TetMesh/small/after_step_vector.step_vector` | [0.857143, 0.857143, 0.857143] | [0.666667, 0.666667, 0.666667] | C2 (resize keeps the extent: new step = extent / cells) | +| `P2TetMesh/small/step_vector` | [0.428571, 0.428571, 0.428571] | [1, 1, 1] | C1 (old array was changed in place by a later setter) | +| `P2TetMesh/small/to_dict.step_vector` | [0.428571, 0.428571, 0.428571] | [1, 1, 1] | C1 (old array was changed in place by a later setter) | +| `P2TetMesh/unit/after_maximum.maximum` | [9, 9, 9] | [8.33333, 8.33333, 8.33333] | C3 (maximum request is enclosed by less than one cell) | +| `P2TetMesh/unit/after_maximum.step_vector` | [1.42857, 1.42857, 1.42857] | [1.33333, 1.33333, 1.33333] | C2 (resize keeps the extent: new step = extent / cells) | +| `P2TetMesh/unit/after_nsteps.step_vector` | [0.714286, 0.714286, 0.714286] | [0.666667, 0.666667, 0.666667] | C2 (resize keeps the extent: new step = extent / cells) | +| `P2TetMesh/unit/after_origin.step_vector` | [1.42857, 1.42857, 1.42857] | [1.33333, 1.33333, 1.33333] | C2 (resize keeps the extent: new step = extent / cells) | +| `P2TetMesh/unit/after_set_nelements.step_vector` | [0.714286, 0.714286, 0.714286] | [0.8, 0.8, 0.8] | C2 (resize keeps the extent: new step = extent / cells) | +| `P2TetMesh/unit/after_step_vector.step_vector` | [1.42857, 1.42857, 1.42857] | [1.33333, 1.33333, 1.33333] | C2 (resize keeps the extent: new step = extent / cells) | +| `P2TetMesh/unit/step_vector` | [0.714286, 0.714286, 0.714286] | [1, 1, 1] | C1 (old array was changed in place by a later setter) | +| `P2TetMesh/unit/to_dict.step_vector` | [0.714286, 0.714286, 0.714286] | [1, 1, 1] | C1 (old array was changed in place by a later setter) | +| `StructuredGrid/aniso/after_maximum.maximum` | [17.4, 10.5455, 11.2857] | [16.5, 8.8, 10] | C3 (maximum request is enclosed by less than one cell) | +| `StructuredGrid/aniso/after_maximum.nsteps` | [7, 9, 7] | [6, 9, 7] | C3 (maximum request is enclosed by less than one cell) | +| `StructuredGrid/aniso/after_maximum.step_vector` | [1.4, 1.81818, 1.71429] | [1.5, 1.6, 1.5] | C2 (resize keeps the extent: new step = extent / cells) | +| `StructuredGrid/aniso/after_nsteps.step_vector` | [0.7, 0.909091, 0.857143] | [0.75, 0.8, 0.75] | C2 (resize keeps the extent: new step = extent / cells) | +| `StructuredGrid/aniso/after_origin.step_vector` | [1.4, 1.81818, 1.71429] | [1.5, 1.6, 1.5] | C2 (resize keeps the extent: new step = extent / cells) | +| `StructuredGrid/aniso/after_set_nelements.step_vector` | [0.7, 0.909091, 0.857143] | [1, 0.888889, 0.9] | C2 (resize keeps the extent: new step = extent / cells) | +| `StructuredGrid/aniso/after_step_vector.step_vector` | [1.4, 1.81818, 1.71429] | [1.5, 1.6, 1.5] | C2 (resize keeps the extent: new step = extent / cells) | +| `StructuredGrid/aniso/step_vector` | [0.7, 0.909091, 0.857143] | [0.5, 2, 1.5] | C1 (old array was changed in place by a later setter) | +| `StructuredGrid/aniso/to_dict.step_vector` | [0.7, 0.909091, 0.857143] | [0.5, 2, 1.5] | C1 (old array was changed in place by a later setter) | +| `StructuredGrid/small/after_maximum.maximum` | [6.71429, 6.71429, 6.71429] | [5.66667, 5.66667, 5.66667] | C3 (maximum request is enclosed by less than one cell) | +| `StructuredGrid/small/after_maximum.nsteps` | [10, 10, 10] | [11, 11, 11] | C3 (maximum request is enclosed by less than one cell) | +| `StructuredGrid/small/after_maximum.step_vector` | [0.857143, 0.857143, 0.857143] | [0.666667, 0.666667, 0.666667] | C2 (resize keeps the extent: new step = extent / cells) | +| `StructuredGrid/small/after_nsteps.step_vector` | [0.428571, 0.428571, 0.428571] | [0.333333, 0.333333, 0.333333] | C2 (resize keeps the extent: new step = extent / cells) | +| `StructuredGrid/small/after_origin.step_vector` | [0.857143, 0.857143, 0.857143] | [0.666667, 0.666667, 0.666667] | C2 (resize keeps the extent: new step = extent / cells) | +| `StructuredGrid/small/after_set_nelements.step_vector` | [0.428571, 0.428571, 0.428571] | [0.4, 0.4, 0.4] | C2 (resize keeps the extent: new step = extent / cells) | +| `StructuredGrid/small/after_step_vector.step_vector` | [0.857143, 0.857143, 0.857143] | [0.666667, 0.666667, 0.666667] | C2 (resize keeps the extent: new step = extent / cells) | +| `StructuredGrid/small/step_vector` | [0.428571, 0.428571, 0.428571] | [1, 1, 1] | C1 (old array was changed in place by a later setter) | +| `StructuredGrid/small/to_dict.step_vector` | [0.428571, 0.428571, 0.428571] | [1, 1, 1] | C1 (old array was changed in place by a later setter) | +| `StructuredGrid/unit/after_maximum.maximum` | [9, 9, 9] | [8.33333, 8.33333, 8.33333] | C3 (maximum request is enclosed by less than one cell) | +| `StructuredGrid/unit/after_maximum.step_vector` | [1.42857, 1.42857, 1.42857] | [1.33333, 1.33333, 1.33333] | C2 (resize keeps the extent: new step = extent / cells) | +| `StructuredGrid/unit/after_nsteps.step_vector` | [0.714286, 0.714286, 0.714286] | [0.666667, 0.666667, 0.666667] | C2 (resize keeps the extent: new step = extent / cells) | +| `StructuredGrid/unit/after_origin.step_vector` | [1.42857, 1.42857, 1.42857] | [1.33333, 1.33333, 1.33333] | C2 (resize keeps the extent: new step = extent / cells) | +| `StructuredGrid/unit/after_set_nelements.step_vector` | [0.714286, 0.714286, 0.714286] | [0.8, 0.8, 0.8] | C2 (resize keeps the extent: new step = extent / cells) | +| `StructuredGrid/unit/after_step_vector.step_vector` | [1.42857, 1.42857, 1.42857] | [1.33333, 1.33333, 1.33333] | C2 (resize keeps the extent: new step = extent / cells) | +| `StructuredGrid/unit/step_vector` | [0.714286, 0.714286, 0.714286] | [1, 1, 1] | C1 (old array was changed in place by a later setter) | +| `StructuredGrid/unit/to_dict.step_vector` | [0.714286, 0.714286, 0.714286] | [1, 1, 1] | C1 (old array was changed in place by a later setter) | +| `TetMesh/aniso/after_maximum.maximum` | [17.4, 10.5455, 11.2857] | [16.5, 8.8, 10] | C3 (maximum request is enclosed by less than one cell) | +| `TetMesh/aniso/after_maximum.nsteps` | [7, 9, 7] | [6, 9, 7] | C3 (maximum request is enclosed by less than one cell) | +| `TetMesh/aniso/after_maximum.step_vector` | [1.4, 1.81818, 1.71429] | [1.5, 1.6, 1.5] | C2 (resize keeps the extent: new step = extent / cells) | +| `TetMesh/aniso/after_nsteps.step_vector` | [0.7, 0.909091, 0.857143] | [0.75, 0.8, 0.75] | C2 (resize keeps the extent: new step = extent / cells) | +| `TetMesh/aniso/after_origin.step_vector` | [1.4, 1.81818, 1.71429] | [1.5, 1.6, 1.5] | C2 (resize keeps the extent: new step = extent / cells) | +| `TetMesh/aniso/after_set_nelements.step_vector` | [0.7, 0.909091, 0.857143] | [1, 0.888889, 0.9] | C2 (resize keeps the extent: new step = extent / cells) | +| `TetMesh/aniso/after_step_vector.step_vector` | [1.4, 1.81818, 1.71429] | [1.5, 1.6, 1.5] | C2 (resize keeps the extent: new step = extent / cells) | +| `TetMesh/aniso/step_vector` | [0.7, 0.909091, 0.857143] | [0.5, 2, 1.5] | C1 (old array was changed in place by a later setter) | +| `TetMesh/aniso/to_dict.step_vector` | [0.7, 0.909091, 0.857143] | [0.5, 2, 1.5] | C1 (old array was changed in place by a later setter) | +| `TetMesh/small/after_maximum.maximum` | [6.71429, 6.71429, 6.71429] | [5.66667, 5.66667, 5.66667] | C3 (maximum request is enclosed by less than one cell) | +| `TetMesh/small/after_maximum.nsteps` | [10, 10, 10] | [11, 11, 11] | C3 (maximum request is enclosed by less than one cell) | +| `TetMesh/small/after_maximum.step_vector` | [0.857143, 0.857143, 0.857143] | [0.666667, 0.666667, 0.666667] | C2 (resize keeps the extent: new step = extent / cells) | +| `TetMesh/small/after_nsteps.step_vector` | [0.428571, 0.428571, 0.428571] | [0.333333, 0.333333, 0.333333] | C2 (resize keeps the extent: new step = extent / cells) | +| `TetMesh/small/after_origin.step_vector` | [0.857143, 0.857143, 0.857143] | [0.666667, 0.666667, 0.666667] | C2 (resize keeps the extent: new step = extent / cells) | +| `TetMesh/small/after_set_nelements.step_vector` | [0.428571, 0.428571, 0.428571] | [0.4, 0.4, 0.4] | C2 (resize keeps the extent: new step = extent / cells) | +| `TetMesh/small/after_step_vector.step_vector` | [0.857143, 0.857143, 0.857143] | [0.666667, 0.666667, 0.666667] | C2 (resize keeps the extent: new step = extent / cells) | +| `TetMesh/small/step_vector` | [0.428571, 0.428571, 0.428571] | [1, 1, 1] | C1 (old array was changed in place by a later setter) | +| `TetMesh/small/to_dict.step_vector` | [0.428571, 0.428571, 0.428571] | [1, 1, 1] | C1 (old array was changed in place by a later setter) | +| `TetMesh/unit/after_maximum.maximum` | [9, 9, 9] | [8.33333, 8.33333, 8.33333] | C3 (maximum request is enclosed by less than one cell) | +| `TetMesh/unit/after_maximum.step_vector` | [1.42857, 1.42857, 1.42857] | [1.33333, 1.33333, 1.33333] | C2 (resize keeps the extent: new step = extent / cells) | +| `TetMesh/unit/after_nsteps.step_vector` | [0.714286, 0.714286, 0.714286] | [0.666667, 0.666667, 0.666667] | C2 (resize keeps the extent: new step = extent / cells) | +| `TetMesh/unit/after_origin.step_vector` | [1.42857, 1.42857, 1.42857] | [1.33333, 1.33333, 1.33333] | C2 (resize keeps the extent: new step = extent / cells) | +| `TetMesh/unit/after_set_nelements.step_vector` | [0.714286, 0.714286, 0.714286] | [0.8, 0.8, 0.8] | C2 (resize keeps the extent: new step = extent / cells) | +| `TetMesh/unit/after_step_vector.step_vector` | [1.42857, 1.42857, 1.42857] | [1.33333, 1.33333, 1.33333] | C2 (resize keeps the extent: new step = extent / cells) | +| `TetMesh/unit/step_vector` | [0.714286, 0.714286, 0.714286] | [1, 1, 1] | C1 (old array was changed in place by a later setter) | +| `TetMesh/unit/to_dict.step_vector` | [0.714286, 0.714286, 0.714286] | [1, 1, 1] | C1 (old array was changed in place by a later setter) | + + +Step 4. Golden update and re-run. `python packages/loop_common/tests/test_structured_support_reference.py --update` wrote 611 arrays. Then `pytest packages/loop_common/tests/test_structured_support_reference.py packages/loop_common/tests/test_structured_grid_2d_reference.py` passed (9 passed). The file `structured_grid_2d_reference.npz` did not change. + +### Numeric survey (S5 resize change) + +Environment: the worktree, `loop_common`, `loop_interpolation` and the visualisation sources, and `pyevtk`, on `PYTHONPATH`. Run after the golden update, one pytest process at a time. + +| Command | Result | Time | +|---|---|---| +| `pytest tests/unit/modelling` | 1540 passed, 2 skipped | 457.78 s | +| `pytest packages/loop_interpolation/tests` | 615 passed, 47 skipped, 2 xfailed | 70.95 s | +| `pytest packages/loop_common/tests` | 406 passed | 69.60 s | +| `pytest tests/unit --ignore=tests/unit/modelling` | 504 passed, 2 skipped | 677.84 s | + +No test failed and no test number changed, so no test file needs a change. The only changed numbers are the 78 golden keys above (file `packages/loop_common/tests/data/structured_support_reference.npz`). Nothing under `packages/loop_interpolation/tests/golden` changed (the 615 tests there pass without a regeneration). No failure is left, and so no root failure from the S4 gate was seen in these four runs. + + +## S5 resize and cell-centre results + +Commands (PowerShell, the worktree `LoopStructural-005-s5` on `PYTHONPATH`): + +- `python -m pytest packages/loop_common/tests/test_resize_paths.py packages/loop_common/tests/test_cell_centre_sampling.py -p no:cacheprovider -q --durations=0` -> 100 passed in 0.19 s. The slowest item took 0.01 s (setup). +- `python -m pytest packages/loop_common/tests/test_resize_paths.py -p no:cacheprovider -q -s -k a_b_a` -> 3 passed (shows the measured drift). +- `python -m pytest packages/loop_common/tests/test_cell_centre_sampling.py -p no:cacheprovider -q` -> 60 passed in 0.12 s. +- `python -m ruff check LoopStructural packages` -> All checks passed. This includes the three `loop_common` source files. + +Counts, as pytest reports them. Resize file: 40 cells (12 + 3 + 3 + 3 + 1 + 1 + 3 + 1 + 10 + 3). The task text said 36. The count differs because `test_2d_counts_are_derived` has 10 cells (2 grids x 5 checks), not 6. Cell-centre file: 60 cells (5 tests x 12). Total 100, which equals the 24 failed + 76 passed of `## S5 pre-fix failures`. + +Comparison with `## S5 pre-fix failures`, test by test. All 24 failed cells now pass. + +| Test function | Before | Now | +|---|---|---| +| `test_nsteps_setter_keeps_extent` | 3 / 3 failed | 3 / 3 pass | +| `test_step_vector_setter_encloses_extent` | 2 / 3 failed (`[G1]`, `[G3]`) | 3 / 3 pass, `[G1]` and `[G3]` now pass | +| `test_nsteps_below_two_is_rejected_without_change` | 1 / 1 failed | 1 / 1 pass | +| `test_nsteps_not_whole_is_rejected_without_change` | 1 / 1 failed | 1 / 1 pass | +| `test_set_nelements_keeps_extent` | 3 / 3 failed | 3 / 3 pass | +| `test_set_nelements_small_request_is_clamped` | 1 / 1 failed | 1 / 1 pass | +| `test_2d_counts_are_derived` | 10 / 10 failed | 10 / 10 pass | +| `test_rectilinear_resize_is_rejected` | 3 / 3 failed | 3 / 3 pass | +| `test_reapplying_same_value_10_times_changes_nothing_3d` | 12 / 12 passed | 12 / 12 pass (coverage) | +| `test_step_vector_a_b_a_extent_drift` | 3 / 3 passed | 3 / 3 pass (coverage) | +| Cell-centre file, 5 tests | 60 / 60 passed | 60 / 60 pass | + +Measured extent drift of the A-B-A test (`step_vector` set to A, then B = 1.5 A, then A). The bound is A + B. + +| Grid | Measured drift | Bound (A + B) | +|---|---|---| +| G1 | [1, 2, 0] | [2.5, 5, 7.5] | +| G2 | [0.5, 0.5, 0.5] | [1.25, 1.25, 1.25] | +| G3 | [10, 100, 7.5] | [25, 250, 18.75] | + +All drifts are inside the bound. The pre-fix drift was G1 [1, 0, 3], G2 [0.5, 0.5, 0.5], G3 [0, 0, 7.5], also inside the bound. The pre-fix and post-fix values differ because the setters now keep the extent and enclose it by less than one cell. + +Evidence for FR-015 and SC-011. `test_centres_are_cell_midpoints_in_world` is the evidence. Each cell of the exported grid (`structured_grid()`) has one centre, as the mean of its 8 corner nodes, and the centres are all different. `test_local_coordinates_grid_associates_own_centres` is a containment check: each cell centre falls inside the `structured_grid(local_coordinates=True)` grid (no -1 from `cell_index_at`). That grid is the axis-aligned enclosure, so a one-to-one count against it does not hold, and the test does not claim it. `test_centres_round_trip_and_index_own_cells` is a consistency check on the cell grid of the box itself, built from the same linspace as the centres. It is not evidence for FR-015 or SC-011. `test_centres_match_exported_cells` and `test_centres_inside_box` are coverage. + +OP12: not reproduced. All 60 cell-centre cells passed before the fix and pass now. `_bounding_box.py` was not changed. + +The named test `test_each_cell_has_exactly_one_centre` of the task text was replaced, because its premise was wrong (see `## S5 pre-fix failures`, note on the one-to-one test). `test_centres_are_cell_midpoints_in_world` is a guard test, because OP12 is not reproduced: it passed before the fix and passes after the fix. + +## S5 diagnostics pre-fix failures + +Method. The test file ran on the unfixed code, in one call: `python -m pytest tests/unit/modelling/test_diagnostics_frames.py -p no:cacheprovider -q -rf --no-header` (PowerShell, the worktree `LoopStructural-005-s5`). No misfit fix commit exists at this point. The result is 147 failed and 0 passed in 142.77 s. No test passes. The file has 147 cells. + +Test file (`test_diagnostics_frames.py`): 147 failed, 0 passed. + +| Test function | Failed / total | First error | +|---|---|---| +| `test_misfit_leaves_state_unchanged` | 24 / 24 | `AssertionError: data['gradient'] changed`. The old methods leave a different array in `data["gradient"]` after the call, in all 24 cells (3 rotations x 4 anchors x 2 engines) | +| `test_new_misfit_methods_leaves_state_unchanged` | 24 / 24 | `AttributeError`: `GeologicalFeature` has no attribute `evaluate_gradient_cosines` | +| `test_gradient_cosine_matches_analytic_direction_exact_field` | 24 / 24 | `assert np.all(np.abs(1.0 - old) <= 1e-10)` is False. The old method gives shape (2, 4) from one column and local vectors. See the values below | +| `test_gradient_cosine_matches_analytic_direction_exact_field_solved_model` | 24 / 24 | `AttributeError`: no attribute `evaluate_gradient_cosines` | +| `test_value_residuals_are_absolute_difference` | 12 / 12 | `AttributeError`: no attribute `evaluate_value_residuals` | +| `test_value_misfit_is_normalised_by_field_range` | 12 / 12 | `AttributeError`: no attribute `evaluate_value_residuals` | +| `test_value_misfit_with_zero_range_does_not_divide` | 2 / 2 | `[equal]`: `RuntimeWarning: divide by zero encountered in divide`. `[nan]`: `assert not np.any(np.isnan(result))` fails, the result has NaN | +| `test_misfit_with_no_rows_is_empty` | 1 / 1 | `AttributeError`: no attribute `evaluate_gradient_cosines` | +| `test_old_misfit_shapes_are_kept` | 4 / 4 | `[gradient-only]`, `[normal-only]`, `[equal-families]`: `AttributeError`, no attribute `evaluate_gradient_cosines`. `[different-families]`: `ValueError`, inhomogeneous shape (the old method gives two lists of different length) | +| `test_stale_misfit_is_not_rebuilt` | 4 / 4 | `[FDI-interpolator]`, `[PLI-interpolator]`: `AssertionError: the call rebuilt the interpolator` (`assert True == False`). `[FDI-builder]`, `[PLI-builder]`: `AssertionError: the call rebuilt the builder` | +| `test_stale_misfit_warns_once_and_returns_nan` | 4 / 4 | `AssertionError: []` (`assert 0 == 1`). The old methods give no `UserWarning` | +| `test_stale_verification_report_does_not_rebuild` | 4 / 4 | `AttributeError`: no attribute `verification_report` | +| `test_report_blocks_are_separate` | 2 / 2 | `AttributeError`: no attribute `verification_report` | +| `test_unsolved_model_gives_no_residual_blocks` | 2 / 2 | `AttributeError`: no attribute `verification_report` | +| `test_exclusions_are_counted_and_explained` | 2 / 2 | `AttributeError`: no attribute `verification_report` | +| `test_gradient_outside_points_give_no_false_unattributed_model` | 1 / 1 | `AttributeError`: no attribute `verification_report` | +| `test_gradient_outside_points_with_pli_are_measured` | 1 / 1 | `AttributeError`: no attribute `verification_report` | + +Group counts: 24 + 24 + 24 + 24 + 12 + 12 + 2 + 1 + 4 + 4 + 4 + 4 + 2 + 2 + 2 + 1 + 1 = 147. + +Items of the task text. + +| Item | Result | +|---|---| +| State test with the old methods | Reproduced. `test_misfit_leaves_state_unchanged` fails in 24 of 24 cells. `data["gradient"]` is a different array after `evaluate_gradient_misfit` (the method divides the stored rows in place) | +| Exact-field cosine test | Reproduced. The field is exact (the gradient is `R1_UNIT` in the world frame), so each cosine must be 1. The old method gives shape (2, 4) and not (8,). Both rows are equal, because the slice `grad[:, :3:6]` is the one column 0 and not the 3 vector columns. The vectors are local, not world. Values of the first row, anchor 0, FDI: rotation 0: [0.2001, 0.2632, 0.1671, 0.2877] (largest `|1 - cos|` 0.833). Rotation 30: [0.0232, 0.0700, -0.0781, 0.1532] (1.078). Rotation 90: [-0.3002, -0.3158, -0.4457, -0.1918] (1.446). The bound is 1e-10 | +| Stale-state test | Reproduced. The call rebuilds the model: `up_to_date` goes from `False` to `True`, in 4 of 4 cells (`[FDI-interpolator]`, `[PLI-interpolator]`: the interpolator flag. `[FDI-builder]`, `[PLI-builder]`: the builder flag) | +| Tests of the new names | Reproduced. 88 cells fail only because the name does not exist: `evaluate_gradient_cosines` 52 cells (24 + 24 + 1 + 3), `verification_report` 12 cells (4 + 2 + 2 + 2 + 1 + 1), `evaluate_value_residuals` 24 cells (12 + 12) | +| `test_value_residuals_are_absolute_difference` | Not reproduced as a value defect. It fails with `AttributeError` (the test uses `evaluate_value_residuals`, which does not exist now). This is not a reproduction of a value defect. The old value path is covered by `test_value_misfit_*` | + +Measured value (a short script in the scratchpad, not saved in the repository). It uses `world_foliation_model` of `tests.fixtures.coordinate_frames`, as `_build` does, with `ANCHORS[1]`, rotation 30 and the FDI engine. + +| Measure | Value | +|---|---| +| Largest absolute change of `interpolator.data["gradient"]` after one call of `evaluate_gradient_misfit` | 8.2254 (the array has 4 rows x 7 columns) | +| Largest change in each column (x, y, z, vx, vy, vz, weight) | 2.9211, 7.5854, 8.2254, 0.5485, 0.2712, 0.6789, 0.9139 | + +The change is large in the position columns and in the vector columns, and the weight column also changes. The old method divides the stored rows in place, so the stored data of the interpolator are not the input data after the call. + +## S5 rollback checks + +Date: 2026-11-05. The check compares the S5 tree with the S4 end code. + +Set-up: + +- The throw-away worktree was at commit 1c029d21 (the docs commit that is the S5 start commit). `git diff --stat 1c029d21 c644c8d7 -- LoopStructural packages` prints nothing. No code differs between 1c029d21 and the S4 end commit c644c8d7, so 1c029d21 holds the S4 end code. +- Each tree had its own `PYTHONPATH`. In each tree, `LoopStructural.__file__` and `loop_interpolation.__file__` pointed into that tree. +- The throw-away worktree was removed with `git worktree remove --force` and `git worktree prune`. `git worktree list` shows no throw-away worktree. + +| Check | Result | Note | +|---|---|---| +| (a) `dataclasses.fields` of the diagnostics classes | pass | `ConstraintFamilyDiagnostics` and `RegionCoverageDiagnostics` are equal. `ConstraintDiagnosticsReport` has only 4 fields on S4 (`interpolator_type`, `families`, `region_coverage`, `outside_model_points`). These 4 are equal on S5. S5 adds `physical_residuals` (the fifth field), `assembled_residuals` and `exclusions`, all with defaults. The task text says "first five", but S4 has only four. | +| (b) `to_dict()` of a report with no verification | pass | One `ConstraintFamilyDiagnostics` and one `RegionCoverageDiagnostics`. The two dicts are equal. | +| (c) `from loop_interpolation import *` | pass | The S4 names are a subset of the S5 names. S5 adds `ExclusionRecord`, `ResidualSummary`, `build_exclusion_records` and `direction_cosines`. No S4 name is missing. | +| (d) `GeologicalFeature.evaluate_gradient_misfit` and `evaluate_value_misfit` | pass | The parameters are equal on both trees (`(self)`). Only the return annotation differs: S5 has `-> 'np.ndarray'`, S4 has none. The result shape for one family (gradient data, normal data dropped, fixture `world_foliation_model`) is `(1, 4)` on both trees. The value misfit has shape `(27,)` on both. No warning on either tree. | +| (e) `StructuredGrid` with an unchanged grid | pass | `origin`, `maximum`, `nsteps` and `step_vector` are equal on both trees (origin [0, 1, 2], nsteps [7, 8, 9], step_vector [1, 2, 0.5], maximum [6, 15, 6]). | +| (f) 2D `maximum` assertions | pass | `pytest packages/loop_common/tests/test_2d_discrete_support.py tests/unit/interpolator/test_2d_discrete_support.py -p no:cacheprovider -q` gives 19 passed on S5 and 19 passed on S4. | + +## S5 post-fix results + +Date: 2026-11-05. The orchestrator ran the gates on commit `6b9e6467`. The `PYTHONPATH` held the worktree `LoopStructural-005-s5`, `packages/loop_common/src`, `packages/loop_interpolation/src`, `packages/loopstructural_visualisation/src` and a scratch folder with `pyevtk`. The root run and the package run were separate calls (PowerShell), because the conftest plugins clash. + +| Command | Result | Time | +|---|---|---| +| `ruff check LoopStructural packages` | All checks passed | - | +| `pytest packages/loop_common/tests packages/loop_interpolation/tests` | 1021 passed, 47 skipped, 2 xfailed, 7 warnings | 140.08 s | +| Root `pytest -W default::DeprecationWarning` | 2216 passed, 4 skipped, 12 warnings | 936.67 s (15 min 36 s) | + +No gate failed. The root run has no failure, and so the S4 failure (`test_write_vol_gocad_is_broken_for_real_bounding_box`) is closed. + +The pytest summary has no library warning with the text "is deprecated" that comes from S5. The log tail (`gate_s5.log`) holds only the last 12 lines of the root run, and it has no "is deprecated" line. + +Cells of the new test files. + +`tests/unit/modelling/test_diagnostics_frames.py` has 147 cells: + +| Test function | Cells | +|---|---| +| `test_misfit_leaves_state_unchanged` | 24 | +| `test_new_misfit_methods_leaves_state_unchanged` | 24 | +| `test_gradient_cosine_matches_analytic_direction_exact_field` | 24 | +| `test_gradient_cosine_matches_analytic_direction_exact_field_solved_model` | 24 | +| `test_value_residuals_are_absolute_difference` | 12 | +| `test_value_misfit_is_normalised_by_field_range` | 12 | +| `test_value_misfit_with_zero_range_does_not_divide` | 2 | +| `test_misfit_with_no_rows_is_empty` | 1 | +| `test_old_misfit_shapes_are_kept` | 4 | +| `test_stale_misfit_is_not_rebuilt` | 4 | +| `test_stale_misfit_warns_once_and_returns_nan` | 4 | +| `test_stale_verification_report_does_not_rebuild` | 4 | +| `test_report_blocks_are_separate` | 2 | +| `test_unsolved_model_gives_no_residual_blocks` | 2 | +| `test_exclusions_are_counted_and_explained` | 2 | +| `test_gradient_outside_points_give_no_false_unattributed_model` | 1 | +| `test_gradient_outside_points_with_pli_are_measured` | 1 | + +The sum is 147. All 147 failed before the fix (`## S5 diagnostics pre-fix failures`) and all pass now. + +`packages/loop_interpolation/tests/test_residual_exclusion_report.py` has 9 tests (no parametrised cells): `test_direction_cosines_pure_function`, `test_report_without_verification_is_unchanged`, `test_builder_records_keep_their_own_key`, `test_unknown_reason_keeps_its_count`, `test_unattributed_rows_are_reported`, `test_gradient_outside_points_give_no_false_unattributed`, `test_dropped_rows_none_gives_no_unattributed`, `test_surfe_style_dropped_equals_outside`, `test_residual_summary_is_importable`. All 9 pass. + +Measured `TOL_SOLVED` = 9e-4. The measured maximum `|1 - cos|` is 4.414e-4 on 24 cells. The margin is a factor of 2. + +Resize and cell-centre files: 100 cells, all pass (`## S5 resize and cell-centre results`). The 3D structured golden file changed in 78 of 611 keys (`## S5 golden review (C-004)`). The 2D golden file did not change. The earlier survey of T078 (modelling 1540 passed, interpolation 615, loop_common 406, tests/unit rest 504) is in the golden section. + +## S5 backend matrix + +| Item | Level | Evidence | +|---|---|---| +| FDI misfit (gradient cosines, value residuals) | mandatory | `test_diagnostics_frames.py` (FDI cells) | +| PLI misfit (gradient cosines, value residuals) | mandatory | `test_diagnostics_frames.py` (PLI cells) | +| Residual and exclusion report, FDI and PLI | mandatory | `test_diagnostics_frames.py`, `test_residual_exclusion_report.py` | +| Resize of the structured grid | mandatory | `test_resize_paths.py` | +| Resize of the 2D grid | mandatory | `test_resize_paths.py` (`test_2d_counts_are_derived`) | +| Resize of the rectilinear grid | mandatory | `test_rectilinear_resize_is_rejected` (resize is rejected with `NotImplementedError`) | +| P2 | not in S5 | S2 and S6 cover it | +| Native engines | not in S5 | S2 and S6 cover them | +| Surfe | not in S5 | S2 and S6 cover it | +| Residual reports for fold and P2 cells | not in S5 | S2 and S6 cover them | + +PLI fact. PLI assembles no gradient rows. In the case with 3 outside points, its gradient family reports an `unattributed` record of 8 (`test_gradient_outside_points_with_pli_are_measured`). The FDI case gives no false `unattributed` record (`test_gradient_outside_points_give_no_false_unattributed_model`). + +## S5 test run time (OP2) + +| Item | Time | +|---|---| +| Root suite (orchestrator, commit `6b9e6467`) | 936.67 s (15 min 36 s). The S4 end run took 805.50 s. | +| Package suites (orchestrator) | 140.08 s. The S4 end run took 150.24 s. | + +The root time is 131 s more than the S4 time. The root suite has 150 more passed tests than at the S4 end (2216 and 2066, and the one S4 failure now passes). The new files add 147 cells to the root suite, and the diagnostics cells are slow (142.77 s on the unfixed code). The two runs had different machine load, so the times do not compare exactly. + +## S5 open points + +- OP12 (cell-centre result): not reproduced. All 60 cell-centre cells passed before and after the fix. The premise of the one-to-one test was corrected. No change to `_bounding_box.py`. +- OP16 (default shape): the default shape is kept. Two families of different length now give the shape (1, N1 + N2). +- OP17 (parity files under `packages/loop_interpolation/tests/golden`): not changed. `git diff --stat 1c029d21 HEAD -- packages/loop_interpolation/tests/golden` prints nothing. diff --git a/specs/README.md b/specs/README.md index b8d31424..2052f714 100644 --- a/specs/README.md +++ b/specs/README.md @@ -9,7 +9,7 @@ not certify that an implementation or a numerical acceptance gate passed. | [002-rbf-interpolator](002-rbf-interpolator/spec.md) | Draft, not reviewed | Continue its existing review and clarification workflow. This cascade did not change it. | | [003-fault-evaluation-performance](003-fault-evaluation-performance/spec.md) | Implementation in progress; 005 contract cascade applied | Its active plan and tasks remain in the existing LoopStructural-native worktree. OP3 and OP4 remain open; completed tasks, including concurrent T043, and measured failures are preserved. | | [004-matrix-free-directional-constraints](004-matrix-free-directional-constraints/spec.md) | Draft; 005 contract cascade applied; task graph checked | Review the updated plan and tasks. OP1 still requires the quadrature measurements before a default is chosen. | -| [005-coordinate-frame-correctness](005-coordinate-frame-correctness/spec.md) | S1 implemented 2026-11-04 on `feat/005-coordinate-frame-correctness` (17 of 17 tasks; evidence in [verification.md](005-coordinate-frame-correctness/verification.md)); S5 and S6 not planned | Review and merge the S1 branch. S2 planned and tasked 2026-11-04 (T018 to T031, waves 8 to 13); run spec-orchestrator for S2. S3 planned and tasked 2026-11-05 ("Stage S3 design" in plan.md; OP4 and OP5 closed; T033 to T050, waves 14 to 22; build after S2 merges). S3 is built on `feat/005-s3-frame-persistence` (pull request #397, stacked on #396; 18 of 18 tasks, evidence in verification.md). S3 merged (#397); a follow-up keeps the legacy load backward compatible (pull request #400). S4 planned and tasked 2026-11-05 ("Stage S4 design" in plan.md; OP8 to OP10 have defaults and block nothing; T051 to T068, waves 23 to 31); run spec-orchestrator for S4. OP1: plan S5 and S6 in later planner runs. OP3: constant-norm accuracy check. OP2 closed 2026-11-04: the R2/R5 matrix stays in the default pytest run. 7 minor review findings are open (run-state.json). | +| [005-coordinate-frame-correctness](005-coordinate-frame-correctness/spec.md) | S1 implemented 2026-11-04 on `feat/005-coordinate-frame-correctness` (17 of 17 tasks; evidence in [verification.md](005-coordinate-frame-correctness/verification.md)); S5 planned 2026-11-05 and revised after review; S6 not planned | Review and merge the S1 branch. S2 planned and tasked 2026-11-04 (T018 to T031, waves 8 to 13); run spec-orchestrator for S2. S3 planned and tasked 2026-11-05 ("Stage S3 design" in plan.md; OP4 and OP5 closed; T033 to T050, waves 14 to 22; build after S2 merges). S3 is built on `feat/005-s3-frame-persistence` (pull request #397, stacked on #396; 18 of 18 tasks, evidence in verification.md). S3 merged (#397); a follow-up keeps the legacy load backward compatible (pull request #400). S4 planned and tasked 2026-11-05 ("Stage S4 design" in plan.md; OP8 to OP10 have defaults and block nothing; T051 to T068, waves 23 to 31); run spec-orchestrator for S4. S5 planned 2026-11-05 and revised the same day after review and five maintainer answers ("Stage S5 design" in plan.md; OP13 to OP15 closed; OP12, OP16 and OP17 have defaults; OP17 blocks the S5 merge only if a parity file changes); S5 tasked 2026-11-05 (T069 to T088, waves 32 to 42; no open point blocks a task; build after S4 merges). OP1: plan S6 in a later planner run. OP3: constant-norm accuracy check. OP2 closed 2026-11-04: the R2/R5 matrix stays in the default pytest run. 7 minor review findings are open (run-state.json). | ## Cascade of 005 diff --git a/tests/unit/modelling/test_diagnostics_frames.py b/tests/unit/modelling/test_diagnostics_frames.py new file mode 100644 index 00000000..80a644f4 --- /dev/null +++ b/tests/unit/modelling/test_diagnostics_frames.py @@ -0,0 +1,535 @@ +"""Frame checks for the misfit methods of a geological feature (U4.4; FR-023; SC-009). + +The data of the foliation are world coordinates. The interpolator stores them in +its local frame. The misfit methods must compare world vectors with world +vectors, and they must not change any state of the model. +""" + +import copy +import itertools +import warnings + +import numpy as np +import pandas as pd +import pytest +from loop_interpolation._diagnostics import UNKNOWN_EXPLANATION, summarise_residuals + +from tests.fixtures.coordinate_frames import ( + ANCHORS, + R1_GRADIENT, + ROTATIONS_DEG, + AnalyticReference, + model_with_frame, + r3_gradient, + r3_value, + world_foliation_model, +) + +ENGINES = ("FDI", "PLI") +R1_UNIT = R1_GRADIENT / np.linalg.norm(R1_GRADIENT) +# Bound of |1 - cosine| for the real solved model. The independent measurement +# over the 24 cells gave a maximum of 4.41e-4. The bound is 2 times this value +# (8.8e-4), rounded up to one significant digit. +TOL_SOLVED = 9e-4 + + +def _cell_id(cell): + anchor, rotation, engine = cell + index = int(np.flatnonzero(np.all(ANCHORS == anchor, axis=1))[0]) + return f"a{index}-r{int(rotation)}-{engine}" + + +CELLS = [ + (ANCHORS[i], rotation, engine) + for i, rotation, engine in itertools.product(range(len(ANCHORS)), ROTATIONS_DEG, ENGINES) +] +CELL_IDS = [_cell_id(cell) for cell in CELLS] +VALUE_CELLS = [ + (ANCHORS[i], rotation, "FDI") + for i, rotation in itertools.product(range(len(ANCHORS)), ROTATIONS_DEG) +] +VALUE_IDS = [_cell_id(cell) for cell in VALUE_CELLS] + +# Row ranges of the data frame of the fixture. +VALUE_ROWS = slice(0, 27) +GRADIENT_ROWS = slice(27, 31) +NORMAL_ROWS = slice(31, 35) + + +def _build(anchor, rotation, engine="FDI", select=None): + """Return a solved model, its foliation and the world data frame.""" + model, feature, data = world_foliation_model(anchor, rotation) + if engine != "FDI" or select is not None: + if select is not None: + data = select(data) + model = model_with_frame(anchor, anchor + (8, 10, 12), rotation) + model.data = data + feature = model.create_and_add_foliation("r1", interpolatortype=engine, nelements=1000) + feature.builder.up_to_date() + assert feature.interpolator.up_to_date + return model, feature, data + + +def _snapshot(feature): + """Return deep copies of the state that a misfit method must not change.""" + interpolator = feature.interpolator + return { + "data": copy.deepcopy(interpolator.data), + "c": copy.deepcopy(interpolator.c), + "up_to_date": interpolator.up_to_date, + "builder_up_to_date": feature.builder._up_to_date, + "weights": { + name: copy.deepcopy(constraint["w"]) + for name, constraint in interpolator.constraints.items() + }, + "report_id": id(interpolator.latest_diagnostics_report), + } + + +def _assert_same_state(before, after): + assert before["data"].keys() == after["data"].keys() + for key, array in before["data"].items(): + assert np.array_equal(array, after["data"][key]), f"data[{key!r}] changed" + assert np.array_equal(before["c"], after["c"]), "c changed" + assert before["up_to_date"] == after["up_to_date"] + assert before["builder_up_to_date"] == after["builder_up_to_date"] + assert before["weights"].keys() == after["weights"].keys() + for name, weights in before["weights"].items(): + assert np.array_equal(weights, after["weights"][name]), f"weights of {name} changed" + assert before["report_id"] == after["report_id"] + + +def _independent_cosines(feature): + """Return the world cosines of the stored rows, with no change to the state. + + The function copies the stored rows and uses only the public frame methods. + """ + interpolator = feature.interpolator + cosines = [] + for key in ("gradient", "normal"): + rows = np.array(interpolator.data[key], dtype=float, copy=True) + if rows.shape[0] == 0: + continue + positions = interpolator.local_to_world_points(rows[:, :3]) + vectors = interpolator.local_to_world_vectors(rows[:, 3:6]) + model_vectors = feature.evaluate_gradient(positions) + dots = np.einsum("ij,ij->i", model_vectors, vectors) + dots /= np.linalg.norm(model_vectors, axis=1) * np.linalg.norm(vectors, axis=1) + cosines.append(dots) + return np.concatenate(cosines) + + +def _world_value_rows(feature, data): + """Return the world points and values of the value rows.""" + rows = data.iloc[VALUE_ROWS] + return rows[["X", "Y", "Z"]].to_numpy(float), rows["val"].to_numpy(float) + + +@pytest.mark.parametrize("cell", CELLS, ids=CELL_IDS) +def test_misfit_leaves_state_unchanged(cell): + """The old misfit methods change no array, flag, weight or report.""" + anchor, rotation, engine = cell + _, feature, _ = _build(anchor, rotation, engine) + before = _snapshot(feature) + feature.evaluate_gradient_misfit() + feature.evaluate_value_misfit() + _assert_same_state(before, _snapshot(feature)) + + +@pytest.mark.parametrize("cell", CELLS, ids=CELL_IDS) +def test_new_misfit_methods_leaves_state_unchanged(cell): + """The new misfit methods and the report change no array, flag, weight or report.""" + anchor, rotation, engine = cell + _, feature, _ = _build(anchor, rotation, engine) + before = _snapshot(feature) + feature.evaluate_gradient_cosines() + feature.evaluate_value_residuals() + feature.verification_report() + _assert_same_state(before, _snapshot(feature)) + + +@pytest.mark.parametrize("cell", CELLS, ids=CELL_IDS) +def test_gradient_cosine_matches_analytic_direction_exact_field(cell, monkeypatch): + """With an exact field the world cosine of each stored row is 1.""" + anchor, rotation, engine = cell + _, feature, _ = _build(anchor, rotation, engine) + + def analytic_gradient(points, *args, **kwargs): + return np.tile(R1_UNIT, (len(np.atleast_2d(points)), 1)) + + monkeypatch.setattr(feature, "evaluate_gradient", analytic_gradient) + old = np.asarray(feature.evaluate_gradient_misfit(), dtype=float) + assert old.shape == (2, 4) + assert np.all(np.abs(1.0 - old) <= 1e-10) + new = feature.evaluate_gradient_cosines() + assert new.shape == (8,) + assert np.all(np.abs(1.0 - new) <= 1e-10) + + +@pytest.mark.parametrize("cell", CELLS, ids=CELL_IDS) +def test_gradient_cosine_matches_analytic_direction_exact_field_solved_model(cell): + """The real solved field has a world cosine near 1 at each stored row.""" + anchor, rotation, engine = cell + _, feature, _ = _build(anchor, rotation, engine) + independent = _independent_cosines(feature) + print(f"max |1 - cos| independent {_cell_id(cell)}: {np.max(np.abs(1.0 - independent)):.3e}") + new = feature.evaluate_gradient_cosines() + assert new.shape == (8,) + assert np.max(np.abs(1.0 - new)) <= TOL_SOLVED + + +@pytest.mark.parametrize("cell", VALUE_CELLS, ids=VALUE_IDS) +def test_value_residuals_are_absolute_difference(cell): + """The residual is the absolute difference at the 27 value rows.""" + anchor, rotation, engine = cell + _, feature, data = _build(anchor, rotation, engine) + points, values = _world_value_rows(feature, data) + expected = np.abs(values - feature.evaluate_value(points)) + residuals = feature.evaluate_value_residuals() + assert residuals.shape == (27,) + assert not np.any(np.isnan(residuals)) + np.testing.assert_allclose(residuals, expected, rtol=0, atol=1e-8) + + +@pytest.mark.parametrize("cell", VALUE_CELLS, ids=VALUE_IDS) +def test_value_misfit_is_normalised_by_field_range(cell): + """The old value misfit is the residual divided by the field range.""" + anchor, rotation, engine = cell + _, feature, _ = _build(anchor, rotation, engine) + expected = feature.evaluate_value_residuals() / (feature.max() - feature.min()) + np.testing.assert_allclose(feature.evaluate_value_misfit(), expected, rtol=0, atol=1e-12) + + +@pytest.mark.parametrize("mode", ["equal", "nan"]) +def test_value_misfit_with_zero_range_does_not_divide(mode, monkeypatch): + """A zero or NaN field range gives the unnormalised difference, with no warning.""" + _, feature, data = _build(ANCHORS[1], 30.0, "FDI") + points, values = _world_value_rows(feature, data) + expected = np.abs(values - feature.evaluate_value(points)) + low, high = (2.0, 2.0) if mode == "equal" else (np.nan, np.nan) + monkeypatch.setattr(feature, "min", lambda: low) + monkeypatch.setattr(feature, "max", lambda: high) + with warnings.catch_warnings(): + warnings.simplefilter("error") + result = feature.evaluate_value_misfit() + assert not np.any(np.isnan(result)) + np.testing.assert_allclose(result, expected, rtol=0, atol=1e-8) + + +def test_misfit_with_no_rows_is_empty(): + """A model with no gradient rows or no value rows gives empty arrays.""" + _, feature, _ = _build(ANCHORS[1], 30.0, "FDI", select=lambda d: d.iloc[VALUE_ROWS]) + cosines = feature.evaluate_gradient_cosines() + assert cosines.shape == (0,) + _, feature, _ = _build( + ANCHORS[1], 30.0, "FDI", select=lambda d: d.iloc[27:35].reset_index(drop=True) + ) + residuals = feature.evaluate_value_residuals() + assert residuals.shape == (0,) + + +def _rows(gradient, normal): + """Return a selector for the value rows and the first rows of each family.""" + + def select(data): + parts = [data.iloc[VALUE_ROWS]] + parts.append(data.iloc[GRADIENT_ROWS].iloc[:gradient]) + parts.append(data.iloc[NORMAL_ROWS].iloc[:normal]) + return pd.concat(parts, ignore_index=True) + + return select + + +@pytest.mark.parametrize( + ("gradient", "normal", "shape"), + [(4, 0, (1, 4)), (0, 4, (1, 4)), (4, 4, (2, 4)), (4, 2, (1, 6))], + ids=["gradient-only", "normal-only", "equal-families", "different-families"], +) +def test_old_misfit_shapes_are_kept(gradient, normal, shape): + """The old method keeps its shapes and holds the values of the new method.""" + _, feature, _ = _build(ANCHORS[1], 30.0, "FDI", select=_rows(gradient, normal)) + old = np.asarray(feature.evaluate_gradient_misfit(), dtype=float) + assert old.shape == shape + np.testing.assert_array_equal(old.ravel(), feature.evaluate_gradient_cosines()) + + +# --------------------------------------------------------------------------- +# Stale model, report blocks and exclusions (U4.4, U4.5; FR-023, FR-024; SC-009) +# --------------------------------------------------------------------------- + +STALE_KINDS = ("interpolator", "builder") +STALE_CELLS = list(itertools.product(ENGINES, STALE_KINDS)) +STALE_IDS = [f"{engine}-{kind}" for engine, kind in STALE_CELLS] +ANCHOR = ANCHORS[1] +ROTATION = 30.0 + + +def _make_stale(feature, kind): + """Make the interpolator or the builder stale, and check the flag. + + The function sets ``builder._up_to_date`` directly. The method + ``set_not_up_to_date`` also makes the interpolator stale. + """ + if kind == "interpolator": + feature.interpolator.up_to_date = False + assert feature.interpolator.up_to_date is False + assert feature.builder._up_to_date is True + else: + feature.builder._up_to_date = False + assert feature.builder._up_to_date is False + assert feature.interpolator.up_to_date is True + + +def _record_exclusion(feature): + """Add one builder record of excluded rows, to make sure that the report has one.""" + feature.builder.excluded_rows["r1_equality"] = {"invalid_reference_value": 2} + + +def _assert_has_recorded_exclusion(report): + """Check that the report holds the record of ``_record_exclusion``.""" + records = [r for r in report.exclusions if r.family == "r1_equality"] + assert len(records) == 1 + assert records[0].reason == "invalid_reference_value" + assert records[0].count == 2 + + +def _stale_feature(engine, kind): + """Return a solved feature that is then made stale.""" + _, feature, _ = _build(ANCHOR, ROTATION, engine) + _record_exclusion(feature) + _make_stale(feature, kind) + return feature + + +def _call_with_warnings(method): + """Call a method and return its result and the UserWarning records.""" + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + result = method() + return result, [w for w in caught if issubclass(w.category, UserWarning)] + + +@pytest.mark.parametrize(("engine", "kind"), STALE_CELLS, ids=STALE_IDS) +def test_stale_misfit_is_not_rebuilt(engine, kind): + """The old misfit methods do not rebuild a stale model (FR-031).""" + feature = _stale_feature(engine, kind) + before = _snapshot(feature) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + feature.evaluate_gradient_misfit() + feature.evaluate_value_misfit() + after = _snapshot(feature) + assert after["up_to_date"] == before["up_to_date"], "the call rebuilt the interpolator" + assert ( + after["builder_up_to_date"] == before["builder_up_to_date"] + ), "the call rebuilt the builder" + _assert_same_state(before, after) + + +@pytest.mark.parametrize(("engine", "kind"), STALE_CELLS, ids=STALE_IDS) +def test_stale_misfit_warns_once_and_returns_nan(engine, kind): + """A stale model gives one UserWarning for each method and NaN values.""" + feature = _stale_feature(engine, kind) + for method, shape in ( + (feature.evaluate_gradient_misfit, (2, 4)), + (feature.evaluate_value_misfit, (27,)), + ): + result, caught = _call_with_warnings(method) + assert len(caught) == 1, [str(w.message) for w in caught] + assert kind in str(caught[0].message).lower() + result = np.asarray(result, dtype=float) + assert result.shape == shape + assert np.all(np.isnan(result)) + for method, count in ( + (feature.evaluate_gradient_cosines, 8), + (feature.evaluate_value_residuals, 27), + ): + result, caught = _call_with_warnings(method) + assert len(caught) == 1, [str(w.message) for w in caught] + assert kind in str(caught[0].message).lower() + assert result.shape == (count,) + assert np.all(np.isnan(result)) + + +@pytest.mark.parametrize(("engine", "kind"), STALE_CELLS, ids=STALE_IDS) +def test_stale_verification_report_does_not_rebuild(engine, kind): + """A stale model gives empty residual blocks, exclusions and one warning.""" + feature = _stale_feature(engine, kind) + before = _snapshot(feature) + report, caught = _call_with_warnings(feature.verification_report) + _assert_same_state(before, _snapshot(feature)) + assert len(caught) == 1, [str(w.message) for w in caught] + assert kind in str(caught[0].message).lower() + assert report.physical_residuals == {} + assert report.assembled_residuals == {} + _assert_has_recorded_exclusion(report) + + +@pytest.mark.parametrize("engine", ENGINES) +def test_report_blocks_are_separate(engine): + """The physical and the assembled residuals are two separate blocks.""" + _, feature, _ = _build(ANCHOR, ROTATION, engine) + report = feature.verification_report() + assert report.physical_residuals is not report.assembled_residuals + assert len(report.physical_residuals) > 0 + assert len(report.assembled_residuals) > 0 + for name, entry in report.physical_residuals.items(): + assert entry.kind == "physical", name + for name, entry in report.assembled_residuals.items(): + assert entry.kind == "assembled", name + + value = report.physical_residuals["value"] + assert value.unit == "field" + expected_rms = float(np.sqrt(np.mean(feature.evaluate_value_residuals() ** 2))) + assert abs(value.rms - expected_rms) <= 1e-12 + + assembled = feature.interpolator.calculate_residual_for_constraints() + assert set(report.assembled_residuals) == set(assembled) + for name, residuals in assembled.items(): + assert report.assembled_residuals[name] == summarise_residuals( + name, "assembled", "assembled", residuals + ) + + payload = report.to_dict() + assert set(payload["physical_residuals"]) == set(report.physical_residuals) + assert set(payload["assembled_residuals"]) == set(report.assembled_residuals) + assert payload["physical_residuals"] is not payload["assembled_residuals"] + assert all(v["kind"] == "physical" for v in payload["physical_residuals"].values()) + assert all(v["kind"] == "assembled" for v in payload["assembled_residuals"].values()) + + +@pytest.mark.parametrize("engine", ENGINES) +def test_unsolved_model_gives_no_residual_blocks(engine): + """A model with no solution has no residual blocks, but it has exclusions.""" + _, feature, _ = _build(ANCHOR, ROTATION, engine) + _record_exclusion(feature) + feature.interpolator.c = None + report = feature.verification_report() + assert report.physical_residuals == {} + assert report.assembled_residuals == {} + _assert_has_recorded_exclusion(report) + + +@pytest.mark.parametrize("engine", ENGINES) +def test_exclusions_are_counted_and_explained(engine): + """A restricted reference gives counted and explained exclusion records.""" + _, feature, _ = _build(ANCHOR, ROTATION, engine) + builder = feature.builder + model = feature.model + reference = AnalyticReference( + lambda p: r3_value(p, ANCHOR), lambda p: r3_gradient(p, ANCHOR), "restricted" + ) + lower, upper = ANCHOR + (2.1, 2.1, 2.1), ANCHOR + (5.9, 7.9, 9.9) + reference.support_box = (lower, upper) + builder.add_equality_constraints(reference, lambda p: np.ones(len(p), dtype=bool)) + builder.add_orthogonal_feature(reference, B=0.0, sample_rule=None) + builder.up_to_date() + assert feature.interpolator.up_to_date + + support = feature.interpolator.support + bounding_box = model.bounding_box + nodes = bounding_box.reproject(support.nodes) + selected = np.all((nodes >= lower) & (nodes <= upper), axis=1) + samples = bounding_box.reproject(support.barycentre) + valid = np.all((samples >= lower) & (samples <= upper), axis=1) + expected = { + "restricted_equality": ("invalid_reference_value", int((~selected).sum())), + "restricted_orthogonal": ("invalid_reference_direction", int((~valid).sum())), + } + report = feature.verification_report() + for family, (reason, count) in expected.items(): + assert count > 0 + records = [r for r in report.exclusions if (r.family, r.reason) == (family, reason)] + assert len(records) == 1, family + assert records[0].count == count + assert records[0].explanation + assert records[0].explanation != UNKNOWN_EXPLANATION + + +def _with_outside_gradient(data): + """Return the data with 3 more gradient rows far outside the model box.""" + extra = data.iloc[GRADIENT_ROWS].iloc[:3].copy() + extra[["X", "Y", "Z"]] = extra[["X", "Y", "Z"]].to_numpy(float) + 100.0 + return pd.concat([data.iloc[VALUE_ROWS], data.iloc[GRADIENT_ROWS], extra], ignore_index=True) + + +def test_gradient_outside_points_give_no_false_unattributed_model(): + """Three gradient points outside the box give outside_model 3, not unattributed.""" + _, feature, _ = _build(ANCHOR, ROTATION, "FDI", select=_with_outside_gradient) + assert feature.interpolator.data["gradient"].shape[0] == 7 + records = feature.verification_report().exclusions + outside = [r for r in records if (r.family, r.reason) == ("gradient", "outside_model")] + assert len(outside) == 1 + assert outside[0].count == 3 + assert not [r for r in records if (r.family, r.reason) == ("gradient", "unattributed")] + + +def test_gradient_outside_points_with_pli_are_measured(): + """For PLI, record the measured counts. PLI assembles no gradient rows.""" + _, feature, _ = _build(ANCHOR, ROTATION, "PLI", select=_with_outside_gradient) + records = feature.verification_report().exclusions + outside = sum(r.count for r in records if (r.family, r.reason) == ("gradient", "outside_model")) + unattributed = sum( + r.count for r in records if (r.family, r.reason) == ("gradient", "unattributed") + ) + print(f"PLI gradient: unattributed={unattributed} outside_model={outside}") + assert outside == 3 + + +def _faulted_built_model(): + """Build a model with one fault and one foliation, then update it.""" + from LoopStructural import GeologicalModel + + model = GeologicalModel([0, 0, 0], [1, 1, 1]) + rows = [ + [0.5, 0.5, 0.5, 1, 0, 0, 0, "fault", 0], + [0.5, 0.5, 0.5, 0, 0, 1, 1, "fault", 0], + [0.5, 0.5, 0.5, 0, 1, 0, 2, "fault", 0], + ] + for x in (0.2, 0.8): + for y in (0.2, 0.8): + rows.append([x, y, 0.3, np.nan, np.nan, np.nan, np.nan, "strati", 0.0]) + rows.append([x, y, 0.7, np.nan, np.nan, np.nan, np.nan, "strati", 1.0]) + rows.append([0.3, 0.3, 0.5, 0, 0, 1, np.nan, "strati", np.nan]) + model.data = pd.DataFrame( + rows, columns=["X", "Y", "Z", "nx", "ny", "nz", "coord", "feature_name", "val"] + ) + fault = model.create_and_add_fault("fault", 0.1, nelements=1e3) + strati = model.create_and_add_foliation("strati", nelements=1e3) + model.update() + return model, fault, strati + + +def test_built_faulted_feature_is_not_stale(): + """A built feature with a fault gives finite misfits and no warning.""" + _, _, strati = _faulted_built_model() + assert len(strati.builder.faults) == 1 + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + cosines = strati.evaluate_gradient_cosines() + residuals = strati.evaluate_value_residuals() + gradient_misfit = strati.evaluate_gradient_misfit() + value_misfit = strati.evaluate_value_misfit() + report = strati.verification_report() + messages = [str(w.message) for w in caught if issubclass(w.category, UserWarning)] + assert messages == [], messages + assert np.all(np.isfinite(cosines)) and np.all(np.isfinite(residuals)) + assert np.all(np.isfinite(gradient_misfit)) and np.all(np.isfinite(value_misfit)) + assert report.physical_residuals and report.assembled_residuals + + +def test_stale_fault_component_is_not_rebuilt(): + """A stale component of a fault gives one warning, NaN, and no rebuild.""" + model, fault, strati = _faulted_built_model() + component = fault.builder.builders[0] + component.set_not_up_to_date(model) + flags = (component._up_to_date, component._interpolator.up_to_date) + assert flags[0] is False + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + residuals = strati.evaluate_value_residuals() + messages = [str(w.message) for w in caught if issubclass(w.category, UserWarning)] + assert len(messages) == 1, messages + assert "builder" in messages[0] or "interpolator" in messages[0] + assert np.all(np.isnan(residuals)) + assert (component._up_to_date, component._interpolator.up_to_date) == flags