diff --git a/COMPAT.md b/COMPAT.md index 98834c540..4e5aa48e0 100644 --- a/COMPAT.md +++ b/COMPAT.md @@ -47,6 +47,7 @@ unchanged, but the new path is preferred going forward. | `FaultBuilder.__init__`/`FoldedFeatureBuilder.__init__`/`StructuralFrameBuilder.__init__` `bounding_box` param annotated as `LoopStructural.geometry._bounding_box.BoundingBox` | annotated as `loop_common.geometry._bounding_box.BoundingBox` — `LoopStructural.geometry.BoundingBox` now re-exports `loop_common.geometry.BoundingBox` (ROADMAP 2c-3 closed); accepted argument type is unchanged, only the class's canonical module path | 2026-07-30 | | `BoundingBox(global_origin=..., global_maximum=...)` (local/global split, `origin`/`maximum` pre-shifted to be near-zero) | `BoundingBox(origin=..., maximum=...)` with `origin`/`maximum` always in world coordinates, plus `set_local_transform(local_origin=...)` for the interpolation frame and `project()`/`reproject()` as a proper affine transform. `global_origin`/`global_maximum` constructor args and properties are removed; `GeologicalModel`'s `scale()`/`rescale()` public methods keep their existing signature and behavior. | 2026-07-30 | | `GeologicalModel.from_file(file)` (always loads via `dill`/`pickle`, no opt-out) | `GeologicalModel.from_file(file, allow_pickle=True)` — same default behavior (still unpickles trusted files with no code change required), but `allow_pickle=False` now refuses to unpickle and raises `LoopValueError` instead, since deserialising an untrusted pickle/dill file can execute arbitrary code. A runtime warning is also now logged whenever pickle-based loading is used. For untrusted/JSON-based input, use `GeologicalModel.from_recipe_dict`/`to_recipe_dict` instead. | 2026-07-30 | +| `create_and_add_fault(fault_center=...)` projected the centre into local coordinates; `fault_buffer`/`force_mesh_geometry` were replaced by a default 0.2 buffer; an inferred centre had z = 0; any 3x3 matrix was accepted as a rotation | `fault_center` is a world point; the requested buffer and forced geometry are kept; an inferred centre keeps its elevation; `set_local_transform` and `rotation_xy` raise `RigidRotationError` for a non-rigid matrix. Extra observations added to work around these defects are no longer necessary (spec 005) | 2026-11-04 | 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/core/_model_feature_factory.py b/LoopStructural/modelling/core/_model_feature_factory.py index ac29bb074..13be2d5dd 100644 --- a/LoopStructural/modelling/core/_model_feature_factory.py +++ b/LoopStructural/modelling/core/_model_feature_factory.py @@ -748,8 +748,6 @@ def build_fault( model._add_faults(fault_frame_builder, features=faults) # add data - if fault_center is not None and ~np.isnan(fault_center).any(): - fault_center = model.scale(fault_center, inplace=False) # Keep the supplied fault-axis values unchanged; the previous self-assignment # was only present to satisfy a linter and did not affect behavior. fault_frame_builder.create_data_from_geometry( @@ -768,8 +766,6 @@ def build_fault( fault_dip_anisotropy=fault_dip_anisotropy, fault_pitch=fault_pitch, ) - if "force_mesh_geometry" not in kwargs: - fault_frame_builder.set_mesh_geometry(kwargs.get("fault_buffer", 0.2), 0) if "splay" in kwargs and "splayregion" in kwargs: fault_frame_builder.add_splay(kwargs["splay"], kwargs["splayregion"]) diff --git a/LoopStructural/modelling/features/builders/_fault_builder.py b/LoopStructural/modelling/features/builders/_fault_builder.py index eeafd9606..0279fc43a 100644 --- a/LoopStructural/modelling/features/builders/_fault_builder.py +++ b/LoopStructural/modelling/features/builders/_fault_builder.py @@ -73,6 +73,8 @@ def __init__( self._fault_major_axis = None self._fault_intermediate_axis = None self._fault_centre = None + self._fault_buffer: float = 0.2 + self._force_mesh_geometry: bool = False self._suspend_rebuild = False def update_geometry(self, points): @@ -102,7 +104,7 @@ def update_geometry(self, points): def _on_geometry_change( self, - force_mesh_geometry: bool = False, + force_mesh_geometry: bool | None = None, fault_buffer: float | None = None, fault_trace_anisotropy: float = 1.0, fault_dip_anisotropy: float = 1.0, @@ -116,6 +118,10 @@ def _on_geometry_change( """ if getattr(self, "raw_data", None) is None or self.raw_data.empty: return + if force_mesh_geometry is None: + force_mesh_geometry = self._force_mesh_geometry + if fault_buffer is None: + fault_buffer = self._fault_buffer # Use a copy of the DataFrame so callers don't observe mutations self.create_data_from_geometry( @@ -129,7 +135,7 @@ def _on_geometry_change( w=1.0, points=False, force_mesh_geometry=force_mesh_geometry, - fault_buffer=fault_buffer if fault_buffer is not None else 0.2, + fault_buffer=fault_buffer, fault_trace_anisotropy=fault_trace_anisotropy, fault_dip=fault_dip if fault_dip is not None else getattr(self, "fault_dip", 90), fault_dip_anisotropy=fault_dip_anisotropy, @@ -148,7 +154,7 @@ def fault_normal_vector(self, v): if len(v.shape) != 1 or v.shape[0] != 3: raise ValueError("fault_normal_vector must be a 3 element array") if v.ndim != 1: - v = v[0,:] + v = v[0, :] if v.shape[0] != 3: raise ValueError("fault_normal_vector must be a 3 element array") arr = np.array(v, dtype=float) @@ -289,6 +295,8 @@ def create_data_from_geometry( Pitch angle of the fault. """ # Suspend automatic rebuilds while we populate internal attributes to avoid recursion + self._fault_buffer = fault_buffer + self._force_mesh_geometry = force_mesh_geometry prev_suspend = getattr(self, "_suspend_rebuild", False) self._suspend_rebuild = True try: @@ -299,37 +307,52 @@ def create_data_from_geometry( ].to_numpy() self.fault_dip = fault_dip if fault_normal_vector is None: - gx_mask = np.logical_and(fault_frame_data["coord"] == 0, fault_frame_data["gx"].notna()) - nx_mask = np.logical_and(fault_frame_data["coord"] == 0, fault_frame_data["nx"].notna()) - value_mask = np.logical_and(fault_frame_data["coord"] == 0, fault_frame_data["val"] == 0) + gx_mask = np.logical_and( + fault_frame_data["coord"] == 0, fault_frame_data["gx"].notna() + ) + nx_mask = np.logical_and( + fault_frame_data["coord"] == 0, fault_frame_data["nx"].notna() + ) + value_mask = np.logical_and( + fault_frame_data["coord"] == 0, fault_frame_data["val"] == 0 + ) if not fault_frame_data.loc[gx_mask].empty: - fault_normal_vector = fault_frame_data.loc[ - gx_mask, - ["gx", "gy", "gz"], - ].to_numpy().mean(axis=0) + fault_normal_vector = ( + fault_frame_data.loc[ + gx_mask, + ["gx", "gy", "gz"], + ] + .to_numpy() + .mean(axis=0) + ) elif not fault_frame_data.loc[nx_mask].empty: - fault_normal_vector = fault_frame_data.loc[ - nx_mask, - ["nx", "ny", "nz"], - ].to_numpy().mean(axis=0) + fault_normal_vector = ( + fault_frame_data.loc[ + nx_mask, + ["nx", "ny", "nz"], + ] + .to_numpy() + .mean(axis=0) + ) else: fault_surface_pts = fault_frame_data.loc[ value_mask, ["X", "Y", "Z"], ].to_numpy() - fault_surface_pts = fault_surface_pts[ - ~np.isnan(fault_surface_pts).any(axis=1) - ] + fault_surface_pts = fault_surface_pts[~np.isnan(fault_surface_pts).any(axis=1)] if fault_surface_pts.shape[0] >= 3: pts_3d = fault_surface_pts - fault_surface_pts.mean(axis=0) cov_matrix = pts_3d.T @ pts_3d eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix) # If points span a surface, use the smallest eigenvector as plane normal - if eigenvalues[-1] > 0 and (eigenvalues[1] / eigenvalues[-1]) > plane_line_threshold: + if ( + eigenvalues[-1] > 0 + and (eigenvalues[1] / eigenvalues[-1]) > plane_line_threshold + ): fault_normal_vector = eigenvectors[:, 0] - + if fault_normal_vector is None: # Fall back to line-on-map strike logic pts = fault_trace - fault_trace.mean(axis=0) @@ -340,9 +363,7 @@ def create_data_from_geometry( strike_vector /= np.linalg.norm(strike_vector) fault_normal_vector = np.cross(strike_vector, [0, 0, 1]) - rotation_matrix = rotation( - strike_vector[None, :], np.array([90 - fault_dip]) - ) + rotation_matrix = rotation(strike_vector[None, :], np.array([90 - fault_dip])) fault_normal_vector = np.einsum( "ijk,ik->ij", rotation_matrix, fault_normal_vector[None, :] )[0] @@ -352,35 +373,66 @@ def create_data_from_geometry( if fault_pitch is not None: rotation_matrix = rotation(fault_normal_vector[None, :], np.array([fault_pitch])) - fault_slip_vector = np.einsum("ijk,ik->ij", rotation_matrix, fault_normal_vector[None, :])[0] + fault_slip_vector = np.einsum( + "ijk,ik->ij", rotation_matrix, fault_normal_vector[None, :] + )[0] if fault_slip_vector is None: - if fault_frame_data.loc[ - np.logical_and(fault_frame_data["coord"] == 1, fault_frame_data["nx"].notna())].shape[0]>0: - fault_slip_vector = fault_frame_data.loc[ - np.logical_and(fault_frame_data["coord"] == 1, fault_frame_data["nx"].notna()), - ["nx", "ny", "nz"], - ].to_numpy().mean(axis=0) + if ( + fault_frame_data.loc[ + np.logical_and( + fault_frame_data["coord"] == 1, fault_frame_data["nx"].notna() + ) + ].shape[0] + > 0 + ): + fault_slip_vector = ( + fault_frame_data.loc[ + np.logical_and( + fault_frame_data["coord"] == 1, fault_frame_data["nx"].notna() + ), + ["nx", "ny", "nz"], + ] + .to_numpy() + .mean(axis=0) + ) else: - fault_slip_vector = np.cross(fault_normal_vector, [1., 0., 0.]) + fault_slip_vector = np.cross(fault_normal_vector, [1.0, 0.0, 0.0]) if np.linalg.norm(fault_slip_vector) == 0: - fault_slip_vector = np.cross(fault_normal_vector, [0., 1., 0.]) + fault_slip_vector = np.cross(fault_normal_vector, [0.0, 1.0, 0.0]) fault_slip_vector /= np.linalg.norm(fault_slip_vector) if fault_center is None: - fault_trace = fault_frame_data.loc[ - np.logical_and(fault_frame_data["coord"] == 0, fault_frame_data["val"] == 0), - ["X", "Y"], - ].to_numpy() - fault_center = fault_trace.mean(axis=0) - fault_center = np.array([fault_center[0], fault_center[1], 0.0]) + centre_rows = fault_frame_data.loc[ + np.logical_and(fault_frame_data["coord"] == 0, fault_frame_data["val"] == 0), + ["X", "Y", "Z"], + ].to_numpy(dtype=float) + fault_center = np.nanmean(centre_rows, axis=0) + for axis, axis_name in enumerate(("X", "Y", "Z")): + if not np.isfinite(fault_center[axis]): + logger.warning( + "Inferred fault centre has no finite %s coordinate; using 0.0", + axis_name, + ) + fault_center[axis] = 0.0 if not isinstance(fault_center, np.ndarray): fault_center = np.array(fault_center) if fault_center.shape[0] != 3: raise ValueError("fault_center must be a 3 element array") # Assign to properties (these won't trigger _on_geometry_change because suspend flag is True) - self.fault_normal_vector = fault_normal_vector / np.linalg.norm(fault_normal_vector) - self.fault_slip_vector = fault_slip_vector / np.linalg.norm(fault_slip_vector) + # Copy the input so the array of the caller does not change. + normal = np.array(fault_normal_vector, dtype=float) + normal /= np.linalg.norm(normal) + slip = np.array(fault_slip_vector, dtype=float) + slip_perp = slip - np.dot(slip, normal) * normal + slip_norm = np.linalg.norm(slip_perp) + if slip_norm < 1e-12: + raise ValueError( + "fault slip vector is parallel to the fault normal; supply a slip direction in the fault plane" + ) + fault_slip_vector = slip_perp / slip_norm + self.fault_normal_vector = normal + self.fault_slip_vector = fault_slip_vector self.fault_centre = fault_center if major_axis is None: @@ -395,6 +447,7 @@ def create_data_from_geometry( self.add_data_from_data_frame(fault_frame_data) self.origin = self.model.bounding_box.origin self.maximum = self.model.bounding_box.maximum + self.set_mesh_geometry(fault_buffer, None) return major_axis = np.max(distance) logger.info(f"Fault major axis using map length: {major_axis}") @@ -404,15 +457,18 @@ def create_data_from_geometry( minor_axis = major_axis / 2.0 if intermediate_axis is None: intermediate_axis = major_axis - logger.info(f"Fault intermediate axis not set, using major axis: {intermediate_axis}") + logger.info( + f"Fault intermediate axis not set, using major axis: {intermediate_axis}" + ) self.fault_minor_axis = minor_axis self.fault_major_axis = major_axis self.fault_intermediate_axis = intermediate_axis self.fault_normal_vector /= np.linalg.norm(self.fault_normal_vector) - fault_slip_vector /= np.linalg.norm(fault_slip_vector) + fault_slip_vector = self.fault_slip_vector.copy() # check if slip vector is inside fault plane, if not project onto fault plane # if not np.isclose(normal_vector @ slip_vector, 0): strike_vector = np.cross(self.fault_normal_vector, fault_slip_vector) + strike_vector /= np.linalg.norm(strike_vector) self.fault_strike_vector = strike_vector fault_edges = np.zeros((2, 3)) @@ -538,6 +594,21 @@ def create_data_from_geometry( ] if major_axis is not None: + fault_frame_data.loc[ + len(fault_frame_data), + ["X", "Y", "Z", "feature_name", "nx", "ny", "nz", "val", "coord", "w"], + ] = [ + fault_center[0], + fault_center[1], + fault_center[2], + self.name, + strike_vector[0] / major_axis, + strike_vector[1] / major_axis, + strike_vector[2] / major_axis, + np.nan, + 2, + w, + ] fault_tips[0, :] = fault_center[:3] + strike_vector * 0.5 * major_axis fault_tips[1, :] = fault_center[:3] - strike_vector * 0.5 * major_axis self.update_geometry(fault_tips) @@ -596,8 +667,10 @@ def create_data_from_geometry( ] self.update_geometry(fault_depth) - # TODO need to add data here - fault_slip_vector /= intermediate_axis + # The depth points above use the unit slip vector. + # The gradient row uses the slip vector divided by intermediate_axis. + # This is a local copy, so self.fault_slip_vector stays a unit vector. + fault_slip_vector = fault_slip_vector / intermediate_axis fault_frame_data.loc[ len(fault_frame_data), [ @@ -704,7 +777,9 @@ def add_fault_trace_anisotropy(self, w: float = 1.0): """ if w > 0: if self.fault_normal_vector is None: - raise ValueError("fault_normal_vector must be initialized before adding anisotropy.") + raise ValueError( + "fault_normal_vector must be initialized before adding anisotropy." + ) plane = np.array([0, 0, 1]) strike_vector = ( @@ -731,7 +806,9 @@ def add_fault_dip_anisotropy(self, w: float = 1.0): """ if w > 0: if self.fault_normal_vector is None: - raise ValueError("fault_normal_vector must be initialized before adding anisotropy.") + raise ValueError( + "fault_normal_vector must be initialized before adding anisotropy." + ) plane = np.array([0, 0, 1]) strike_vector = ( diff --git a/LoopStructural/modelling/features/builders/_geological_feature_builder.py b/LoopStructural/modelling/features/builders/_geological_feature_builder.py index 2ed03f535..ebd6fe527 100644 --- a/LoopStructural/modelling/features/builders/_geological_feature_builder.py +++ b/LoopStructural/modelling/features/builders/_geological_feature_builder.py @@ -83,6 +83,9 @@ def __init__( ) self._orthogonal_features = {} self._equality_constraints = {} + # Rows not installed from reference features: + # {constraint name: {reason: number of rows}} + self.excluded_rows: dict[str, dict[str, int]] = {} # add default parameters self.update_build_arguments( { @@ -305,6 +308,8 @@ def add_data_to_interpolator(self, constrained=False, force_constrained=False, * def install_gradient_constraint(self): if issubclass(type(self.interpolator), DiscreteInterpolator): + for key in [k for k in self.excluded_rows if k.endswith("_orthogonal")]: + del self.excluded_rows[key] for g in self._orthogonal_features.values(): feature, w, _region, step, B, sample_rule = g if w == 0: @@ -323,16 +328,35 @@ def install_gradient_constraint(self): support.nodes, support.elements[tets], sample_rule ) weight = w * np.sqrt(q) + name = f"{feature.name}_orthogonal" + # The sample points are local. The reference feature + # works in world coordinates. + world_points = self.interpolator.local_to_world_points(points) # don't worry about masking for unconformities here - vector = feature.evaluate_gradient(points, ignore_regions=True) + vector = np.asarray( + feature.evaluate_gradient(world_points, ignore_regions=True), dtype=float + ) norm = np.linalg.norm(vector, axis=1) - vector[norm > 0] /= norm[norm > 0, None] + valid = np.all(np.isfinite(vector), axis=1) & np.isfinite(norm) & (norm > 0) + n_invalid = int((~valid).sum()) + if n_invalid > 0: + self.excluded_rows.setdefault(name, {})[ + "invalid_reference_direction" + ] = n_invalid + if not valid.any(): + logger.warning( + f"No valid reference gradients for {name} in {self.name}: " + f"{n_invalid} rows excluded" + ) + continue + vector = vector[valid] / norm[valid, None] + local_vector = self.interpolator.world_to_local_vectors(vector) self.interpolator.add_gradient_orthogonal_constraints( - points, - vector, - w=weight, + points[valid], + local_vector, + w=weight[valid] if np.ndim(weight) > 0 else weight, b=B, - name=f'{feature.name}_orthogonal', + name=name, ) def add_equality_constraints(self, feature, region, scalefactor=1.0): @@ -341,20 +365,32 @@ def add_equality_constraints(self, feature, region, scalefactor=1.0): self._up_to_date = False def install_equality_constraints(self): + for key in [k for k in self.excluded_rows if k.endswith("_equality")]: + del self.excluded_rows[key] for e in self._equality_constraints.values(): + name = f"{e[0].name}_equality" + idc = None try: # assume all parts of structural frame have the same support support = self.interpolator.support # work out the values of the nodes where we want hard - # constraints - idc = np.arange(0, support.n_nodes)[e[1](support.nodes)] - val = e[0].evaluate_value(support.nodes[e[1](support.nodes), :]) - mask = ~np.isnan(val) - self.interpolator.add_equality_constraints(idc[mask], val[mask] * e[2]) + # constraints. The nodes are local. The region and the + # reference feature work in world coordinates. + world_nodes = self.interpolator.local_to_world_points(support.nodes) + selected = np.asarray(e[1](world_nodes), dtype=bool) + idc = np.flatnonzero(selected) + val = np.asarray(e[0].evaluate_value(world_nodes[idc]), dtype=float) + mask = np.isfinite(val) + n_invalid = int((~mask).sum()) + if n_invalid > 0: + self.excluded_rows.setdefault(name, {})["invalid_reference_value"] = n_invalid + self.interpolator.add_equality_constraints(idc[mask], val[mask] * e[2], name=name) except (AttributeError, TypeError, ValueError, RuntimeError) as exc: logger.error(f"Could not add equality for {self.name}") logger.error(f"Exception: {exc}") + count = idc.size if idc is not None else self.interpolator.support.n_nodes + self.excluded_rows.setdefault(name, {})["evaluation_error"] = int(count) def get_value_constraints(self): """ @@ -475,13 +511,10 @@ def set_interpolation_geometry(self, origin, maximum, rotation=None): return # origin/maximum are given in world coordinates (e.g. straight from - # model.bounding_box or fault-frame data); project into the - # interpolator's local frame before writing to the support. Exact - # for translation-only transforms -- no code sets a non-identity - # rotation on the bounding box today. + # model.bounding_box or fault-frame data). Use the local box that + # encloses all corners of the world box. if self.interpolator.bounding_box is not None: - origin = self.interpolator.bounding_box.project(origin) - maximum = self.interpolator.bounding_box.project(maximum) + origin, maximum = self.interpolator.bounding_box.local_enclosure(origin, maximum) self.interpolator.support.origin = origin self.interpolator.support.maximum = maximum @@ -499,16 +532,14 @@ def check_interpolation_geometry(self, data, buffer=0.3): """ if self.interpolator.support is not None: logger.info(f"Checking interpolation geometry for {self.name}") - origin = self.interpolator.support.origin - maximum = self.interpolator.support.maximum - pts = self.model.bounding_box.with_buffer(buffer).regular_grid(local=True) + # The grid and the fault restoration are in world coordinates. + pts = self.model.bounding_box.with_buffer(buffer).regular_grid(local=False) for f in self.faults: pts = f.apply_to_points(pts) - - origin[origin > np.min(pts, axis=0)] = np.min(pts, axis=0)[origin > np.min(pts, axis=0)] - maximum[maximum < np.max(pts, axis=0)] = np.max(pts, axis=0)[ - maximum < np.max(pts, axis=0) - ] + pts = self.interpolator.world_to_local_points(pts) + # Use copies. Do not change the support arrays in place. + origin = np.minimum(self.interpolator.support.origin.copy(), np.min(pts, axis=0)) + maximum = np.maximum(self.interpolator.support.maximum.copy(), np.max(pts, axis=0)) self.interpolator.support.origin = origin self.interpolator.support.maximum = maximum self.update_build_arguments({'nelements': self.interpolator.n_elements}) diff --git a/LoopStructural/modelling/features/fault/_fault_segment.py b/LoopStructural/modelling/features/fault/_fault_segment.py index 1bffa2db0..944675d59 100644 --- a/LoopStructural/modelling/features/fault/_fault_segment.py +++ b/LoopStructural/modelling/features/fault/_fault_segment.py @@ -126,15 +126,9 @@ def fault_ellipsoid(self, **kwargs): try: import pyvista as pv - if self.model is None: - pts = self.fault_centre[None, :] - else: - pts = self.model.rescale(self.fault_centre[None, :], inplace=False) - # pts = self.fault_centre[None, :] + # The fault centre is a world point. Do not convert it again. + pts = self.fault_centre[None, :] fault_ellipsoid = pv.PolyData(pts) - # fault_ellipsoid = pv.PolyData( - # self.model.rescale(self.fault_centre[None, :], inplace=False) - # ) fault_ellipsoid["norm"] = self.fault_normal_vector[None, :] fault_ellipsoid['norm'] /= np.linalg.norm(fault_ellipsoid['norm'], axis=1)[:, None] geom = pv.ParametricEllipsoid( diff --git a/packages/loop_common/src/loop_common/geometry/__init__.py b/packages/loop_common/src/loop_common/geometry/__init__.py index 4e47b8308..53f8ee33c 100644 --- a/packages/loop_common/src/loop_common/geometry/__init__.py +++ b/packages/loop_common/src/loop_common/geometry/__init__.py @@ -1,5 +1,6 @@ from ._bounding_box import BoundingBox from ._point import ValuePoints, VectorPoints +from ._rigid import RigidRotationError, validate_rigid_rotation from ._structured_grid import PropertyGrid, StructuredGrid from ._structured_grid_2d import StructuredGrid2DGeometry from ._structured_grid_3d import StructuredGrid3DGeometry @@ -9,6 +10,7 @@ __all__ = [ "BoundingBox", "PropertyGrid", + "RigidRotationError", "StructuredGrid", "StructuredGrid2DGeometry", "StructuredGrid3DGeometry", @@ -17,4 +19,5 @@ "UnstructuredMeshGeometry", "ValuePoints", "VectorPoints", + "validate_rigid_rotation", ] diff --git a/packages/loop_common/src/loop_common/geometry/_bounding_box.py b/packages/loop_common/src/loop_common/geometry/_bounding_box.py index 3b5f5f5a6..5be4bd598 100644 --- a/packages/loop_common/src/loop_common/geometry/_bounding_box.py +++ b/packages/loop_common/src/loop_common/geometry/_bounding_box.py @@ -1,6 +1,7 @@ from __future__ import annotations import copy +from itertools import product import numpy as np @@ -9,6 +10,7 @@ # from LoopStructural.utils.exceptions import LoopValueError from loop_common.math import rng +from ._rigid import validate_rigid_rotation from ._structured_grid import PropertyGrid logger = getLogger(__name__) @@ -18,7 +20,6 @@ class LoopValueError(ValueError): """Custom error for invalid values in LoopStructural.""" - class BoundingBox: def __init__( self, @@ -74,7 +75,9 @@ def _coerce_point(point, name): if nsteps is not None: if len(nsteps) != dimensions: - logger.warning(f"Nsteps has {len(nsteps)} dimensions but bounding box has {dimensions}") + logger.warning( + f"Nsteps has {len(nsteps)} dimensions but bounding box has {dimensions}" + ) raise LoopValueError("Nsteps has incorrect number of dimensions") if np.any(np.asarray(nsteps) <= 0): raise LoopValueError("Nsteps must be positive integers") @@ -146,17 +149,15 @@ def set_local_transform( """ if local_origin is None: local_origin = np.zeros(self.dimensions, dtype=float) - local_origin = np.asarray(local_origin, dtype=float) + local_origin = np.array(local_origin, dtype=float, copy=True) if local_origin.shape != (self.dimensions,): raise LoopValueError("Local origin has incorrect number of dimensions") + if not np.all(np.isfinite(local_origin)): + raise LoopValueError("Local origin must contain only finite values") if rotation_matrix is None: rotation_matrix = np.eye(self.dimensions, dtype=float) - rotation_matrix = np.asarray(rotation_matrix, dtype=float) - if rotation_matrix.shape != (self.dimensions, self.dimensions): - raise LoopValueError( - f"Rotation matrix must have shape ({self.dimensions}, {self.dimensions})" - ) + rotation_matrix = validate_rigid_rotation(rotation_matrix, self.dimensions) self._local_origin = local_origin self._local_rotation = rotation_matrix @@ -168,6 +169,50 @@ def set_local_transform( self._world_to_local = world_to_local self._local_to_world = np.linalg.inv(world_to_local) + def local_enclosure( + self, + origin: np.ndarray | None = None, + maximum: np.ndarray | None = None, + ) -> tuple[np.ndarray, np.ndarray]: + """Return local bounds enclosing a world-space axis-aligned box. + + Parameters + ---------- + origin : np.ndarray, optional + World-space minimum corner. Defaults to the bounding box origin. + maximum : np.ndarray, optional + World-space maximum corner. Defaults to the bounding box maximum. + + Returns + ------- + tuple[np.ndarray, np.ndarray] + Local-space minimum and maximum coordinates. + + Raises + ------ + LoopValueError + If either bound is missing, invalid, or has the wrong shape. + """ + if origin is None: + origin = self._origin + if maximum is None: + maximum = self._maximum + if origin is None or maximum is None: + raise LoopValueError("Origin and maximum must be set") + + origin = np.asarray(origin, dtype=float) + maximum = np.asarray(maximum, dtype=float) + if origin.shape != (self.dimensions,) or maximum.shape != (self.dimensions,): + raise LoopValueError("Bounds have incorrect number of dimensions") + if not np.all(np.isfinite(origin)) or not np.all(np.isfinite(maximum)): + raise LoopValueError("Bounds must contain only finite values") + if np.any(origin > maximum): + raise LoopValueError("Maximum must be greater than or equal to origin") + + corners = np.array(list(product(*zip(origin, maximum))), dtype=float) + local_corners = self.project(corners) + return np.min(local_corners, axis=0), np.max(local_corners, axis=0) + def _apply_affine( self, xyz: np.ndarray, matrix: np.ndarray, inplace: bool = False ) -> np.ndarray: diff --git a/packages/loop_common/src/loop_common/geometry/_rigid.py b/packages/loop_common/src/loop_common/geometry/_rigid.py new file mode 100644 index 000000000..bd186e087 --- /dev/null +++ b/packages/loop_common/src/loop_common/geometry/_rigid.py @@ -0,0 +1,50 @@ +"""Validation helpers for rigid rotations.""" + +from __future__ import annotations + +import numpy as np + + +class RigidRotationError(ValueError): + """Raised when a matrix is not a proper rotation.""" + + +def validate_rigid_rotation(matrix: np.ndarray, dimensions: int, atol: float = 1e-10) -> np.ndarray: + """Validate and return a copy of a proper rotation matrix. + + Parameters + ---------- + matrix : np.ndarray + Candidate rotation matrix. + dimensions : int + Required matrix dimensions. + atol : float, default=1e-10 + Absolute tolerance for orthonormality and determinant checks. + + Returns + ------- + np.ndarray + Floating point copy of the validated rotation matrix. + + Raises + ------ + RigidRotationError + If the matrix is not a finite proper rotation. + """ + rotation = np.array(matrix, dtype=float) + suffix = "; supply a proper rotation (orthonormal, det +1); scale and shear are not supported" + + if rotation.shape != (dimensions, dimensions): + raise RigidRotationError(f"matrix shape must be ({dimensions}, {dimensions}){suffix}") + if not np.all(np.isfinite(rotation)): + raise RigidRotationError(f"matrix entries must be finite{suffix}") + + orthogonality_error = np.max(np.abs(rotation.T @ rotation - np.eye(dimensions))) + if orthogonality_error > atol: + raise RigidRotationError(f"matrix must be orthonormal{suffix}") + + determinant_error = abs(np.linalg.det(rotation) - 1.0) + if determinant_error > atol: + raise RigidRotationError(f"matrix determinant must be +1{suffix}") + + return rotation 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 b78ac033c..cf1c35146 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,11 +1,12 @@ """Pure 3D regular grid geometry: origin/nsteps/step_vector indexing.""" - import numpy as np from loop_common.logging import get_logger as getLogger from loop_common.utils import LoopException +from ._rigid import validate_rigid_rotation + logger = getLogger(__name__) @@ -87,11 +88,41 @@ def rotation_xy(self, rotation_xy): if rotation_xy is None: return if isinstance(rotation_xy, (float, int)): - rotation_xy = np.array([[np.cos(np.deg2rad(rotation_xy)), -np.sin(np.deg2rad(rotation_xy)), 0], [np.sin(np.deg2rad(rotation_xy)), np.cos(np.deg2rad(rotation_xy)), 0], [0, 0, 1]]) - rotation_xy = np.array(rotation_xy) - if rotation_xy.shape != (3, 3): - raise ValueError(f"Rotation matrix should be 3x3, not {rotation_xy.shape}") - self._rotation_xy = rotation_xy + rotation_xy = np.array( + [ + [np.cos(np.deg2rad(rotation_xy)), -np.sin(np.deg2rad(rotation_xy)), 0], + [np.sin(np.deg2rad(rotation_xy)), np.cos(np.deg2rad(rotation_xy)), 0], + [0, 0, 1], + ] + ) + self._rotation_xy = validate_rigid_rotation(rotation_xy, 3) + + def _cells_for_length(self, length: np.ndarray) -> np.ndarray: + """Return cell counts that cover the requested axis lengths. + + Parameters + ---------- + length : np.ndarray + Requested support length along each grid axis. + + Returns + ------- + np.ndarray + Number of cells required along each axis. + + Raises + ------ + ValueError + If any requested cell count is negative. + """ + cells = np.ceil(length / self.step_vector - 1e-9).astype(np.int64) + cells[cells == 0] = 2 + if np.any(cells < 0): + logger.error( + f"Cannot resize the grid. The proposed number of cells is {cells}, these must be all >= 0" + ) + raise ValueError("Cannot resize the grid.") + return cells @property def step_vector(self): @@ -112,13 +143,8 @@ def origin(self): @origin.setter def origin(self, origin): origin = np.array(origin) - length = self.maximum - origin - length /= self.step_vector - self._nsteps = np.ceil(length).astype(np.int64) - self._nsteps[self._nsteps == 0] = 3 - if np.any(~(self._nsteps > 0)): - logger.error(f"Cannot resize the grid. The proposed number of steps is {self._nsteps}, these must be all > 0") - raise ValueError("Cannot resize the grid.") + maximum = self.maximum + self._nsteps = self._cells_for_length(maximum - origin) + 1 self._origin = origin self._on_geometry_change() @@ -129,13 +155,7 @@ def maximum(self): @maximum.setter def maximum(self, maximum): maximum = np.array(maximum, dtype=float) - length = maximum - self.origin - length /= self.step_vector - self._nsteps = np.ceil(length).astype(np.int64) - self._nsteps[self._nsteps == 0] = 3 - if np.any(~(self._nsteps > 0)): - logger.error(f"Cannot resize the grid. The proposed number of steps is {self._nsteps}, these must be all > 0") - raise ValueError("Cannot resize the grid.") + self._nsteps = self._cells_for_length(maximum - self.origin) + 1 self._on_geometry_change() @property @@ -223,16 +243,26 @@ def cell_corner_indexes(self, cell_indexes: np.ndarray) -> np.ndarray: xcorner = np.array([0, 1, 0, 1, 0, 1, 0, 1]) ycorner = np.array([0, 0, 1, 1, 0, 0, 1, 1]) zcorner = np.array([0, 0, 0, 0, 1, 1, 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[:, :, 2] = cell_indexes[:, None, 2] + corner_indexes[:, :, 2] + zcorner[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, :] + ) + corner_indexes[:, :, 2] = ( + cell_indexes[:, None, 2] + corner_indexes[:, :, 2] + zcorner[None, :] + ) return corner_indexes def global_index_to_cell_index(self, global_index): cell_indexes = np.zeros((global_index.shape[0], 3), dtype=np.int64) cell_indexes[:, 0] = global_index % self.nsteps_cells[0, None] - cell_indexes[:, 1] = (global_index // self.nsteps_cells[0, None]) % self.nsteps_cells[1, None] - cell_indexes[:, 2] = global_index // self.nsteps_cells[0, None] // self.nsteps_cells[1, None] + cell_indexes[:, 1] = (global_index // self.nsteps_cells[0, None]) % self.nsteps_cells[ + 1, None + ] + cell_indexes[:, 2] = ( + global_index // self.nsteps_cells[0, None] // self.nsteps_cells[1, None] + ) return cell_indexes def global_index_to_node_index(self, global_index): diff --git a/packages/loop_common/tests/data/structured_support_reference.npz b/packages/loop_common/tests/data/structured_support_reference.npz index 43a4b0795..f4330e176 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_rigid_transform_contract.py b/packages/loop_common/tests/test_rigid_transform_contract.py new file mode 100644 index 000000000..01c3dd0ba --- /dev/null +++ b/packages/loop_common/tests/test_rigid_transform_contract.py @@ -0,0 +1,154 @@ +"""Contract tests for rigid coordinate transforms and grid resize stability.""" + +import numpy as np +import pytest +from loop_common.geometry import ( + BoundingBox, + RigidRotationError, + StructuredGrid3DGeometry, + validate_rigid_rotation, +) + + +@pytest.mark.parametrize("dimensions", [2, 3]) +@pytest.mark.parametrize("angle", [0.0, 30.0, 90.0]) +def test_validate_rigid_rotation_accepts_proper_rotations(dimensions, angle): + """Accept proper rotations in two and three dimensions.""" + radians = np.deg2rad(angle) + rotation = np.eye(dimensions) + rotation[:2, :2] = [[np.cos(radians), -np.sin(radians)], [np.sin(radians), np.cos(radians)]] + + result = validate_rigid_rotation(rotation, dimensions) + + np.testing.assert_allclose(result, rotation) + assert result is not rotation + + +@pytest.mark.parametrize( + "matrix", + [ + np.diag([2.0, 1.0, 1.0]), + np.array([[1.0, 0.2, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]), + np.diag([1.0, 1.0, 0.0]), + np.diag([-1.0, 1.0, 1.0]), + np.array([[np.nan, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]), + np.array([[np.inf, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]), + np.eye(2), + ], + ids=["scale", "shear", "singular", "reflection", "nan", "infinity", "wrong-shape"], +) +def test_validate_rigid_rotation_rejects_invalid(matrix): + """Reject transforms that are not finite proper rotations.""" + assert issubclass(RigidRotationError, ValueError) + with pytest.raises(RigidRotationError, match=""): + validate_rigid_rotation(matrix, 3) + + +@pytest.mark.parametrize( + "anchor", + [ + (0.0, 0.0, 0.0), + (100.0, 200.0, 300.0), + (-100.0, -200.0, -300.0), + (500000.0, 6500000.0, 100.0), + ], +) +@pytest.mark.parametrize("dimensions", [2, 3]) +@pytest.mark.parametrize("angle", [0.0, 30.0, 90.0]) +def test_point_and_direction_round_trips(anchor, dimensions, angle): + """Round-trip points and vectors without translation leaking into vectors.""" + anchor = np.asarray(anchor[:dimensions], dtype=float) + box = BoundingBox( + origin=anchor - np.array([100.0, 200.0, 50.0][:dimensions]), + maximum=anchor + np.array([100.0, 200.0, 50.0][:dimensions]), + dimensions=dimensions, + ) + radians = np.deg2rad(angle) + rotation = np.array([[np.cos(radians), -np.sin(radians)], [np.sin(radians), np.cos(radians)]]) + if dimensions == 3: + rotation = np.pad(rotation, ((0, 1), (0, 1))) + rotation[2, 2] = 1.0 + box.set_local_transform(local_origin=np.array(anchor), rotation_matrix=rotation) + offset = np.array([3.0, -7.0, 11.0][:dimensions]) + points = np.array([anchor, anchor + offset]) + vectors = np.array([[1.0, 2.0, 3.0][:dimensions], [-4.0, 5.0, -6.0][:dimensions]]) + length = np.linalg.norm(box.maximum - box.origin) + world_magnitude = max(np.linalg.norm(box.origin), np.linalg.norm(box.maximum)) + point_tolerance = 1e-10 * length + 32 * np.finfo(float).eps * max(world_magnitude, length) + + np.testing.assert_allclose( + box.reproject(box.project(points)), points, atol=point_tolerance, rtol=0.0 + ) + np.testing.assert_allclose( + box.reproject(box.project(points[0])), points[0], atol=point_tolerance, rtol=0.0 + ) + np.testing.assert_allclose( + box.reproject_vectors(box.project_vectors(vectors)), vectors, atol=1e-12, rtol=0.0 + ) + np.testing.assert_allclose( + box.reproject_vectors(box.project_vectors(vectors[0])), + vectors[0], + atol=1e-12, + rtol=0.0, + ) + + +@pytest.mark.parametrize("angle", [0.0, 30.0, 90.0]) +def test_local_enclosure_contains_projected_corners(angle): + """Enclose every projected world corner within the returned local bounds.""" + origin = np.array([-10.0, 20.0, -30.0]) + maximum = np.array([40.0, 80.0, 60.0]) + box = BoundingBox(origin=origin, maximum=maximum) + radians = np.deg2rad(angle) + rotation = np.array( + [ + [np.cos(radians), -np.sin(radians), 0.0], + [np.sin(radians), np.cos(radians), 0.0], + [0.0, 0.0, 1.0], + ] + ) + box.set_local_transform(local_origin=origin, rotation_matrix=rotation) + local_minimum, local_maximum = box.local_enclosure() + tolerance = 1e-10 * np.linalg.norm(maximum - origin) + 32 * np.finfo(float).eps * np.max( + maximum - origin + ) + + assert np.all(box.project(box.corners) >= local_minimum - tolerance) + assert np.all(box.project(box.corners) <= local_maximum + tolerance) + + +def test_invalid_rotation_keeps_prior_state(): + """Keep the prior transform when a rotation update fails validation.""" + box = BoundingBox(origin=np.zeros(3), maximum=np.ones(3)) + old_transform = box.world_to_local_matrix + with pytest.raises(RigidRotationError): + box.set_local_transform(rotation_matrix=np.diag([2.0, 1.0, 1.0])) + np.testing.assert_array_equal(box.world_to_local_matrix, old_transform) + + grid = StructuredGrid3DGeometry(rotation_xy=30.0) + old_rotation = grid.rotation_xy.copy() + with pytest.raises(RigidRotationError): + grid.rotation_xy = np.diag([-1.0, 1.0, 1.0]) + np.testing.assert_array_equal(grid.rotation_xy, old_rotation) + + +def test_reapply_unchanged_bounds_is_stable(): + """Keep geometry unchanged when setting its current bounds repeatedly.""" + grid = StructuredGrid3DGeometry( + origin=np.array([-10.0, 20.0, 30.0]), + nsteps=np.array([11, 21, 31]), + step_vector=np.array([2.0, 3.0, 4.0]), + ) + origin = grid.origin.copy() + maximum = grid.maximum.copy() + step_vector = grid.step_vector.copy() + nsteps = grid.nsteps.copy() + + for _ in range(10): + grid.origin = origin.copy() + grid.maximum = maximum.copy() + + np.testing.assert_array_equal(grid.origin, origin) + np.testing.assert_array_equal(grid.maximum, maximum) + np.testing.assert_array_equal(grid.step_vector, step_vector) + np.testing.assert_array_equal(grid.nsteps, nsteps) diff --git a/packages/loop_interpolation/src/loop_interpolation/_discrete_interpolator.py b/packages/loop_interpolation/src/loop_interpolation/_discrete_interpolator.py index 4fd073517..6cb366f87 100644 --- a/packages/loop_interpolation/src/loop_interpolation/_discrete_interpolator.py +++ b/packages/loop_interpolation/src/loop_interpolation/_discrete_interpolator.py @@ -418,7 +418,9 @@ def add_constraints_to_least_squares(self, A, B, idc, w=1.0, name="undefined"): def add_gradient_orthogonal_constraints( self, points: np.ndarray, vectors: np.ndarray, w: float = 1.0 ): - pass + """`points` are target-local support points; `vectors` are target-local + unit directions. + """ def get_regularisation_sample_points(self) -> np.ndarray: raise NotImplementedError( @@ -629,6 +631,9 @@ def add_inequality_feature( Evaluate the feature at the interpolation nodes. Can provide a boolean mask to restrict to only some parts + `feature` receives world points and returns one value for + each support node. + Parameters ---------- feature : BaseFeature @@ -642,7 +647,7 @@ def add_inequality_feature( # flag lower determines whether the feature is a lower bound or upper bound # mask is just a boolean array determining which nodes to apply it to - value = feature(self.support.nodes) + value = feature(self.local_to_world_points(self.support.nodes)) if mask is None: mask = np.ones(value.shape[0], dtype=bool) l = np.zeros(value.shape[0]) - np.inf diff --git a/packages/loop_interpolation/src/loop_interpolation/_geological_interpolator.py b/packages/loop_interpolation/src/loop_interpolation/_geological_interpolator.py index 48e19723c..e6f08a0ad 100644 --- a/packages/loop_interpolation/src/loop_interpolation/_geological_interpolator.py +++ b/packages/loop_interpolation/src/loop_interpolation/_geological_interpolator.py @@ -3,6 +3,7 @@ This module contains the abstract base class for all geological interpolators used in LoopStructural geological modelling framework. """ + from __future__ import annotations import json @@ -209,9 +210,7 @@ def check_array(self, array: np.ndarray): except (TypeError, ValueError) as e: raise LoopTypeError(str(e)) - def _coerce_value_constraint( - self, points: np.ndarray | ValueConstraint - ) -> ValueConstraint: + def _coerce_value_constraint(self, points: np.ndarray | ValueConstraint) -> ValueConstraint: if isinstance(points, ValueConstraint): return points return ValueConstraint.from_array(points, dimensions=self.dimensions) @@ -285,6 +284,88 @@ def _project_constraint_array(self, array: np.ndarray, has_vector: bool = False) array[:, d : 2 * d] = self._project_vectors(array[:, d : 2 * d]) return array + def world_to_local_points(self, points: np.ndarray) -> np.ndarray: + """Convert point locations from world to local coordinates. + + Parameters + ---------- + points : np.ndarray + World coordinates, shape ``(N, d)`` or ``(d,)``. + + Returns + ------- + np.ndarray + Local coordinates with the same shape as ``points``. The input is + returned unchanged if no ``bounding_box`` is set. + + Notes + ----- + The conversion applies the translation and the rotation of the + bounding box. + """ + return self._project_points(points) + + def local_to_world_points(self, points: np.ndarray) -> np.ndarray: + """Convert point locations from local to world coordinates. + + Parameters + ---------- + points : np.ndarray + Local coordinates, shape ``(N, d)`` or ``(d,)``. + + Returns + ------- + np.ndarray + World coordinates with the same shape as ``points``. The input is + returned unchanged if no ``bounding_box`` is set. + + Notes + ----- + The conversion applies the rotation and the translation of the + bounding box. + """ + return self._reproject_points(points) + + def world_to_local_vectors(self, vectors: np.ndarray) -> np.ndarray: + """Convert direction vectors from world to local coordinates. + + Parameters + ---------- + vectors : np.ndarray + World direction vectors, shape ``(N, d)`` or ``(d,)``. + + Returns + ------- + np.ndarray + Local direction vectors with the same shape as ``vectors``. The + input is returned unchanged if no ``bounding_box`` is set. + + Notes + ----- + Vectors use the rotation only. The translation does not apply. + """ + return self._project_vectors(vectors) + + def local_to_world_vectors(self, vectors: np.ndarray) -> np.ndarray: + """Convert direction vectors from local to world coordinates. + + Parameters + ---------- + vectors : np.ndarray + Local direction vectors, shape ``(N, d)`` or ``(d,)``. + + Returns + ------- + np.ndarray + World direction vectors with the same shape as ``vectors``. The + input is returned unchanged if no ``bounding_box`` is set. + + Notes + ----- + Vectors use the rotation only. The translation does not apply. + """ + return self._reproject_vectors(vectors) + @abstractmethod def set_region(self, **kwargs): """Set the interpolation region. diff --git a/packages/loop_interpolation/tests/golden/test_golden.py b/packages/loop_interpolation/tests/golden/test_golden.py index 6dd5e22e7..7ceac8c33 100644 --- a/packages/loop_interpolation/tests/golden/test_golden.py +++ b/packages/loop_interpolation/tests/golden/test_golden.py @@ -136,17 +136,21 @@ def test_matches_golden_field(case, engine, request): if golden is None: pytest.fail("no golden file: run with --update-golden and review the figures first") # The solver stops at a relative residual of 1e-8. Other platforms (macOS - # arm64) add in a different order and differ by about 1e-6, so 1e-6 is - # too tight. A real change of the field is much larger than 1e-5. + # arm64) add in a different order and differ by about 1e-6 in the value, so + # 1e-6 is too tight. The gradient is a derivative of the field, so it + # amplifies that noise by about 1/h (h = 1/8 on the P2 mesh): it needs a + # looser limit than the value. A real change of the field is much larger + # than 1e-4. + dv = np.nanmax(np.abs(r["value"] - golden["value"])) + dg = np.nanmax(np.abs(r["gradient"] - golden["gradient"])) ok = np.allclose( r["value"], golden["value"], atol=1e-5, rtol=0, equal_nan=True - ) and np.allclose(r["gradient"], golden["gradient"], atol=1e-5, rtol=0, equal_nan=True) + ) and np.allclose(r["gradient"], golden["gradient"], atol=1e-4, rtol=0, equal_nan=True) if not ok: from . import render # needs matplotlib fig = render.render_diff(case, engine, r["points"], golden, r) - dv = np.nanmax(np.abs(r["value"] - golden["value"])) - pytest.fail(f"field changed, max |dv| = {dv:.3g}. Diff figure: {fig}") + pytest.fail(f"field changed, max |dv| = {dv:.3g}, max |dgrad| = {dg:.3g}. Diff figure: {fig}") def _compare_family(name, now, ref): diff --git a/specs/005-coordinate-frame-correctness/plan.md b/specs/005-coordinate-frame-correctness/plan.md new file mode 100644 index 000000000..04a2cdd9f --- /dev/null +++ b/specs/005-coordinate-frame-correctness/plan.md @@ -0,0 +1,661 @@ +# Implementation plan: Coordinate and fault-frame correctness + +**Feature**: `005-coordinate-frame-correctness` | **Date**: 2026-11-04 | **Spec**: [spec.md](spec.md) + + + +## Open points + +- **OP1**: This plan gives the exact design only for stage S1 (shared + transform, reference-constraint boundary and fault frame). The spec-planner + limit is 10 code files, and the maintainer selected a staged plan on + 2026-11-04. Stages S2 to 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, after S1 merges. Blocks: the + tasks of S2 to S6. The spec-tasks skill must make tasks for S1 only. +- **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 + requested cells). R5 needs 3 buffers x 3 rotations x 2 centre modes x 2 + engines = 36 support builds. SC-018 forbids skips of mandatory cells. Owner: + the maintainer. Decision needed: run the full matrix in the default + `pytest` run, or put it behind a pytest marker that a named CI job runs. + Default until an answer: the full matrix runs in the default `pytest` run. + The S1 test task measures the time and records it in `verification.md`. Blocks: + nothing in S1 design; it can change the test layout only. + +## Answers + +- 2026-11-04 (wave 3, user decision): the T006 resize repair changes the + golden file `packages/loop_common/tests/data/structured_support_reference.npz`. + 27 values change: `after_origin.nsteps` and `after_maximum.nsteps` and + `after_maximum.maximum`, for 3 support classes and 3 grids. The old values + did not enclose the requested bounds (FR-014 defect). The file is added to + the `Files` of T006. T006 regenerates it with + `python packages/loop_common/tests/test_structured_support_reference.py --update`. + T015 records the change under `## Golden review (C-004)`. +- 2026-11-04 (spec review, user decision): `FaultSegment.fault_ellipsoid` + called `model.rescale` on `builder.fault_centre`, which S1 makes a world + point. `LoopStructural/modelling/features/fault/_fault_segment.py` is added + to the `Files` of T013. The ellipsoid uses the world centre with no + `rescale`. A regression test goes in `test_fault_frame_analytic.py` (T011). +- 2026-11-04 (spec review, user decision): fix two minor record findings in + S1. The golden file keeps only the 27 FR-014 changes; the 33 other keys + get their old values back. `verification.md` gets the correct FR ranges in + its "Pre-fix failures" titles. + +- 2026-11-04 (OP2, maintainer decision, spec C-007): the R2 and R5 matrix + stays in the default `pytest` run with no marker. Measured S1 time: about + 5 minutes for the four new test files. OP2 is closed. + +## Summary + +Observations, fault geometry, masks and public queries stay in world +coordinates. Local assembly and supports use the computational frame of the +target interpolator (`interpolator.bounding_box`). Each boundary converts +points once and directions once (rotation only). + +S1 repairs the shared boundary and the fault frame: + +1. A rigid-rotation check and a corner-enclosure helper in common geometry. +2. Support resize that keeps the node and cell counts consistent. +3. Public world/local conversion methods on `GeologicalInterpolator`. +4. Reference sampling (orthogonality, equality, feature inequality, restored + support expansion) at world locations, with directions converted into the + target basis. +5. A fault factory that does not project an explicit centre, keeps the + requested buffer and forced geometry, infers a 3D centre, and adds one + extent-direction row at the centre. + +S2 to S6 repair fold paths, persistence and lifecycle, display consumers, +grid and diagnostics, and the auxiliary utilities and validation matrix +(see "Stages and rollback"). + +## Technical context + +- **Packages changed**: `LoopStructural` (modelling builders and the feature + factory), `packages/loop_common` (geometry), `packages/loop_interpolation` + (base and discrete interpolators). S1 does not change + `packages/loop_interpolation_core`. +- **Languages and versions**: Python only. Ruff target py39. New modules use + `from __future__ import annotations` for `X | None` hints. +- **New dependencies**: none. +- **Public API and import changes**: additive only in S1. New public names: + `loop_common.geometry.validate_rigid_rotation`, + `BoundingBox.local_enclosure`, + `GeologicalInterpolator.world_to_local_points`, + `GeologicalInterpolator.local_to_world_points`, + `GeologicalInterpolator.world_to_local_vectors`, + `GeologicalInterpolator.local_to_world_vectors`, + `GeologicalFeatureBuilder.excluded_rows`. The signatures of + `GeologicalModel.create_and_add_fault`, `FaultBuilder` and + `set_interpolation_geometry` do not change. Behaviour changes with known + coordinate meaning (FR-028): an explicit `fault_center` is used as a world + point; `fault_buffer` and `force_mesh_geometry` keep their effect; an + inferred centre keeps its elevation; `BoundingBox.set_local_transform` and + `rotation_xy` reject a non-rigid matrix. +- **Saved-data format changes**: none in S1. S3 adds a format version (FR-016). +- **Coordinate contract**: + + | Boundary | Input space | Output space | Function | + |---|---|---|---| + | Public evaluation | world points | world values, world directions | `GeologicalInterpolator.evaluate_value`, `evaluate_gradient` (unchanged) | + | Data to interpolator | world constraint arrays | local stored arrays | `GeologicalInterpolator.set_*_constraints` (unchanged) | + | Support samples to reference feature | local support points | world points | `GeologicalInterpolator.local_to_world_points` | + | Reference direction to assembly | world direction | target-local direction | `GeologicalInterpolator.world_to_local_vectors` | + | Region and mask callbacks | world points | boolean | caller callbacks (documented contract) | + | `DiscreteInterpolator.add_gradient_orthogonal_constraints` | local points, local unit vectors | rows | internal local contract, now documented | + | `DiscreteInterpolator.add_inequality_feature` feature callback | world points | values | documented feature-facing contract | + | Support geometry request | world box | local enclosure | `BoundingBox.local_enclosure` | + | Fault restoration for support expansion | world points | world points | `FaultSegment.apply_to_points` (unchanged) | + + The geological-frame values (fault coordinates 0, 1, 2) are scalars. No + translation applies to them. +- **Performance target**: none specified. OP2 records the test run time. +- **Platforms and backends**: S1 covers Python FDI and PLI/P1. P2, fold, + constant-norm, native and Surfe cells belong to S2 and S6. The native + adapter rejects rotated supports (`_native.py:170`); S1 does not change it. + +## Constitution check + +| Project rule | Result | How the plan agrees | +|---|---|---| +| Work stays within the existing library and package responsibilities | Pass | Changes stay in `loop_common` geometry, `loop_interpolation` base classes and `LoopStructural` builders. One new private module `_rigid.py` in `loop_common.geometry`. No new package. | +| Development commands come from .dev-tools/project.json and the constitution | Pass | Gates use `pytest` and `ruff check LoopStructural packages`, and `black .` for format. | +| Both pytest and ruff check LoopStructural packages must pass before merge | Pass | Each stage row in "Stages and rollback" names both gates. | +| Python changes use numpy-style docstrings and type hints; Ruff targets py39 | Pass | Each new function has a numpy docstring and hints; `from __future__ import annotations` keeps py39 syntax. | +| Black uses line length 100 and preserves string quoting | Pass | Run `black .`; no string normalisation. | +| ROADMAP.md owns architecture stages and decisions; update it when either changes | Pass | S1 adds no ROADMAP stage. S3 records the FR-019 legacy-load decision in ROADMAP.md. The final stage adds a "005 complete" line. | +| Stages can be tested and reversed independently | Pass | Six stages, each with an acceptance check and a rollback target. | +| The public GeologicalModel API stays compatible; a break needs a ROADMAP.md decision and a major-version note | Pass | No signature or return-shape change. Corrected behaviour of known-meaning inputs gets a `COMPAT.md` "Migration notices" entry. S3 handles the FR-019 exception. | +| The QGIS plugin stays in its separate repository | Pass | No plugin file changes. | +| Feature documents stay under specs/NNN-name/ | Pass | `plan.md` and `verification.md` are in `specs/005-coordinate-frame-correctness/`. | +| Releases use release-please with per-package paths | Pass | Conventional commits per package path; no release config change. | +| A constitution amendment requires a pull request and maintainer acceptance | Not applicable | No amendment. | + +## Design + +S1 has the exact design below. S2 to S6 have scope only (OP1). + +### Types and data + +1. **`FaultBuilder` state** (`LoopStructural/modelling/features/builders/_fault_builder.py`). + New private fields, set in `__init__`: + - `self._fault_buffer: float = 0.2` — last requested support buffer, as a + fraction of the longest box side. + - `self._force_mesh_geometry: bool = False` — last forced-geometry request. + + `create_data_from_geometry` writes both fields from its arguments before + any other work. `_on_geometry_change` reads them when its arguments are + `None`. + +2. **`GeologicalFeatureBuilder.excluded_rows`** + (`LoopStructural/modelling/features/builders/_geological_feature_builder.py`). + Type `dict[str, dict[str, int]]`. Key 1: constraint name + (`f"{feature.name}_orthogonal"` or `f"{feature.name}_equality"`). Key 2: + reason. Value: number of rows not installed. Reasons: + - `"invalid_reference_value"` — the reference value is not finite. + - `"invalid_reference_direction"` — the reference gradient is not finite + or has zero length. + - `"evaluation_error"` — the reference evaluation raised an exception; the + count is the number of selected rows. + + The builder resets the dictionary at the start of + `install_gradient_constraint` (orthogonal keys) and + `install_equality_constraints` (equality keys). S5 connects it to the + diagnostics report. + +3. **Generated fault rows** (data frame columns `X, Y, Z, feature_name, nx, + ny, nz, val, coord, w`, world coordinates). S1 adds one row when + `major_axis` is not `None`: + + | X, Y, Z | nx, ny, nz | val | coord | w | + |---|---|---|---|---| + | fault centre (world) | `strike_unit / major_axis` | NaN | 2 | `w` | + + `strike_unit = cross(n, s_perp) / |cross(n, s_perp)|`, where `n` is the + unit normal and `s_perp` is the unit slip direction after the normal + component is removed. For R2, `n=(1,0,0)`, `s=(0,0,-1)`, so + `strike_unit=(0,1,0)` and the row is `(0, 1/20, 0)`. This is the only + automatic direction row (R4, SC-013). Existing tip rows (`val=+-0.5`, + `coord=2`) do not change. + +### Functions and public API + +All names below were checked with `git grep -n -w ` on 2026-11-04 in +`LoopStructural`, `packages` and `tests`: 0 matches for each new name. + +**`packages/loop_common/src/loop_common/geometry/_rigid.py`** (new) + +```python +def validate_rigid_rotation( + matrix: np.ndarray, dimensions: int, atol: float = 1e-10 +) -> np.ndarray: + """Return a float copy of ``matrix`` if it is a proper rigid rotation. + + Raises RigidRotationError (a ValueError subclass in this module) when: + shape != (dimensions, dimensions); any entry is not finite; + max|R^T R - I| > atol; |det(R) - 1| > atol. + The message names the failed condition and the action: + "supply a proper rotation (orthonormal, det +1); scale and shear are + not supported". + """ +``` + +`RigidRotationError(ValueError)` is defined in the same module. It is a +`ValueError` subclass, so existing `except ValueError` callers still work. +Export both names from `loop_common.geometry.__init__`. The module imports +only `numpy`, so it causes no import cycle. + +**`BoundingBox`** (`packages/loop_common/src/loop_common/geometry/_bounding_box.py`) + +- `set_local_transform(local_origin=None, rotation_matrix=None)`: call + `validate_rigid_rotation(rotation_matrix, self.dimensions)` before any field + changes. Also reject a non-finite `local_origin` with `LoopValueError`. + On error the old transform stays (SC-012). +- New method: + + ```python + def local_enclosure( + self, + origin: np.ndarray | None = None, + maximum: np.ndarray | None = None, + ) -> tuple[np.ndarray, np.ndarray]: + """Return (local_min, local_max) that enclose the projection of all + 2**dimensions corners of the world box [origin, maximum]. + Defaults: self.origin and self.maximum. Raises LoopValueError if + either bound is None or not finite, or if origin > maximum on any axis. + """ + ``` + + Implementation: build the corners with `itertools.product` over the two + bounds per axis, call `self.project(corners)`, return `min` and `max` on + axis 0. + +**`StructuredGrid3DGeometry`** (`packages/loop_common/src/loop_common/geometry/_structured_grid_3d.py`) + +- `rotation_xy` setter: after the existing float/int conversion, call + `validate_rigid_rotation(rotation_xy, 3)`. `None` keeps the current value. +- New private method used by both setters: + + ```python + def _cells_for_length(self, length: np.ndarray) -> np.ndarray: + """Number of cells per axis that cover ``length`` with the current step. + + cells = ceil(length / step_vector - 1e-9). An axis with 0 cells gets 2 + cells (the existing degenerate rule of 3 nodes). Raises ValueError + (existing message "Cannot resize the grid.") if any cell count < 0. + """ + ``` + +- `origin` setter: `self._nsteps = self._cells_for_length(self.maximum - origin) + 1`. + Compute `self.maximum` before the origin changes. +- `maximum` setter: + `self._nsteps = self._cells_for_length(maximum - self.origin) + 1`. + + This repairs the 1-cell loss (node count was set to the cell count). + Reapplying unchanged bounds keeps origin, maximum, step and counts + (SC-011, first part). + +**`GeologicalInterpolator`** (`packages/loop_interpolation/src/loop_interpolation/_geological_interpolator.py`) + +New public methods. Each is a documented wrapper of the existing private +method. The private methods stay for compatibility. + +```python +def world_to_local_points(self, points: np.ndarray) -> np.ndarray: ... +def local_to_world_points(self, points: np.ndarray) -> np.ndarray: ... +def world_to_local_vectors(self, vectors: np.ndarray) -> np.ndarray: ... +def local_to_world_vectors(self, vectors: np.ndarray) -> np.ndarray: ... +``` + +Input shape `(N, d)` or `(d,)`; output has the same shape. Vectors use +rotation only. With `bounding_box is None` each method returns its input +(existing behaviour). + +**`DiscreteInterpolator`** (`packages/loop_interpolation/src/loop_interpolation/_discrete_interpolator.py`) + +- `add_gradient_orthogonal_constraints` (abstract): docstring states the + internal contract — `points` are target-local, `vectors` are target-local + unit directions. No signature change. +- `add_inequality_feature(feature, lower=True, mask=None)`: call + `feature(self.local_to_world_points(self.support.nodes))`. Docstring: + `feature` receives world points and returns one value for each node. No + signature change. + +**`GeologicalFeatureBuilder`** (`LoopStructural/modelling/features/builders/_geological_feature_builder.py`) + +- `install_gradient_constraint()`: + 1. `local_pts` = barycentres or quadrature points (existing code). + 2. `world_pts = self.interpolator.local_to_world_points(local_pts)`. + 3. `vector = feature.evaluate_gradient(world_pts, ignore_regions=True)`. + 4. `valid = all(isfinite(vector), axis=1) & (norm > 0)`. + Record `(~valid).sum()` under `"invalid_reference_direction"`. + 5. Normalise `vector[valid]`, then + `local_vec = self.interpolator.world_to_local_vectors(vector[valid])`. + 6. Call `add_gradient_orthogonal_constraints(local_pts[valid], local_vec, + w=weight[valid] if weight is an array else weight, b=B, name=...)`. + 7. If `valid.sum() == 0`, install no rows and log a warning with the count. +- `install_equality_constraints()`: + 1. `world_nodes = self.interpolator.local_to_world_points(support.nodes)`. + 2. `selected = np.asarray(region(world_nodes), dtype=bool)`; + `idc = np.flatnonzero(selected)`. + 3. `val = feature.evaluate_value(world_nodes[idc])`; `valid = isfinite(val)`. + 4. Record `(~valid).sum()` under `"invalid_reference_value"`. + 5. `self.interpolator.add_equality_constraints(idc[valid], + val[valid] * scalefactor, name=f"{feature.name}_equality")`. + 6. The existing `except` clause stays. It also records `idc.size` (or the + node count if `idc` is not set) under `"evaluation_error"`. +- `set_interpolation_geometry(origin, maximum, rotation=None)`: replace the + two-point projection with + `origin, maximum = self.interpolator.bounding_box.local_enclosure(origin, maximum)` + when `bounding_box` is not `None`. The NaN checks and the `dof < 100` loop + do not change. +- `check_interpolation_geometry(data, buffer=0.3)`: + 1. `pts = self.model.bounding_box.with_buffer(buffer).regular_grid(local=False)` + (world). + 2. For each fault: `pts = f.apply_to_points(pts)` (world in, world out). + 3. `local = self.interpolator.world_to_local_points(pts)`. + 4. `origin = np.minimum(support.origin.copy(), local.min(axis=0))`; + `maximum = np.maximum(support.maximum.copy(), local.max(axis=0))`. + 5. Assign `support.origin = origin`, then `support.maximum = maximum`. + The copies stop the in-place change of the support arrays. + +**`FaultBuilder`** (`LoopStructural/modelling/features/builders/_fault_builder.py`) + +- `create_data_from_geometry(...)`: no signature change. + - Store `self._fault_buffer = fault_buffer` and + `self._force_mesh_geometry = force_mesh_geometry` first. + - Inferred centre (`fault_center is None`): use the rows with `coord == 0` + and `val == 0`, columns `X, Y, Z`; `fault_center = np.nanmean(pts, axis=0)`. + If a component is not finite, set it to `0.0` and log a warning that names + the axis. This keeps the elevation (FR-010). + - Slip direction: after normalisation, `s_perp = s - (s @ n) * n`. If + `|s_perp| < 1e-12`, raise `ValueError("fault slip vector is parallel to + the fault normal; supply a slip direction in the fault plane")`. Use + `s_perp / |s_perp|` for `self.fault_slip_vector`, for the strike vector, + and for the depth rows. + - Add the extent-direction row (see "Types and data", item 3). + - Early-return branch (no fault trace): call + `self.set_mesh_geometry(fault_buffer, None)` before `return`. This keeps + the support that the factory previously set. +- `_on_geometry_change(force_mesh_geometry=None, fault_buffer=None, ...)`: + the default of `force_mesh_geometry` changes from `False` to `None`. `None` + means "use the stored value". This is a private method. + +**`create_and_add_fault`** (`LoopStructural/modelling/core/_model_feature_factory.py`) + +- Delete the call `fault_center = model.scale(fault_center, inplace=False)`. + An explicit centre stays a world point (FR-010, U1.1). +- Delete the block + `if "force_mesh_geometry" not in kwargs: fault_frame_builder.set_mesh_geometry(kwargs.get("fault_buffer", 0.2), 0)`. + `create_data_from_geometry` already applies the named `fault_buffer` and + `force_mesh_geometry` (FR-012, U1.4, U1.5). + +### Messages and saved data + +None in S1. S3 owns the saved-format change (FR-016 to FR-019). + +### Behaviour + +Call sequence for a fault through the public API after S1: + +1. `GeologicalModel.create_and_add_fault(...)` calls the factory with world + data and an optional world `fault_center`. +2. The factory makes `FaultBuilder` with `model.bounding_box` and calls + `create_data_from_geometry` with the unchanged centre. +3. The builder makes the world rows, sets `origin`/`maximum` (forced or from + rows), and calls `set_mesh_geometry(self._fault_buffer, None)`. +4. `set_mesh_geometry` calls `set_interpolation_geometry` on the three + coordinate builders. Each builder projects the world box with + `local_enclosure` and resizes its support with the repaired setters. +5. `build` sends world rows to `set_*_constraints`. The interpolator converts + them once into local arrays (existing code). +6. `install_gradient_constraint` samples the reference coordinate at world + points, converts the directions into the target basis, and installs rows + at local points. +7. Evaluation and displacement take world points and return world values and + directions (existing code). + +Errors: an invalid rotation raises `RigidRotationError` before any state +change. A slip parallel to the normal raises `ValueError` before the data +frame is added to the builder. + +Cache invalidation: S1 changes no cache. S3 owns transform changes after +population (FR-020). + +### Stages and rollback + +| Stage | Files and behaviour | Acceptance check | Rollback target and compatibility check | +|---|---|---|---| +| S1 | Shared transform boundary, reference relationships for FDI and PLI, fault frame, buffer and forced geometry, rotation check, resize stability (this design) | New S1 tests pass for FDI and PLI; `pytest`; `ruff check LoopStructural packages`; `pytest packages/loop_common/tests packages/loop_interpolation/tests` | Revert the S1 commits to master `ebc7c699`. Check: the existing `tests/unit/modelling/test__fault_builder.py` and `packages/loop_common/tests/test_bounding_box.py` pass on the reverted tree. No saved format changed. | +| 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) | SC-007, SC-008, SC-013 for every mandatory fold and P2 cell | Revert to the S1 merge. Check: S1 tests still pass. | +| S3 | Persistence format version, transform and coordinate-space metadata; 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) | 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. | +| 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 + +Shared rule for every stage (FR-031): write each regression test first, run +it on the unfixed tree, and record the failure (test id, observed value) in +`specs/005-coordinate-frame-correctness/verification.md`. Then fix and record +the pass. A finding of the 1.7.0 investigation that does not reproduce on the +current tree is recorded as "not reproduced"; its test stays as a guard. + +- **Unit and regression tests (S1)**: + - `packages/loop_common/tests/test_rigid_transform_contract.py`: + `validate_rigid_rotation` accepts rotations of 0, 30 and 90 degrees about + z; rejects scale, shear, singular, reflection (det -1), NaN/inf and wrong + shape (2D and 3D). SC-001 point round trips (`Tp`) and direction round + trips (`1e-12`) for 1 point and arrays, for all 4 anchors. `local_enclosure` + contains all 8 projected corners for each rotation. Invalid rotation: + `set_local_transform` and `rotation_xy` leave the prior state unchanged + (SC-012, rotation part). Reapply unchanged support bounds 10 times: 0 + changes to origin, maximum, step and counts (SC-011, resize part). + - `tests/unit/modelling/test_reference_constraint_frames.py`: R3 for FDI + and PLI, 4 anchors x 3 rotations. Equality region `x' > 4` selects 168 of + 336 nodes unrotated; for rotations the expected set comes from the world + nodes. 0 missing or duplicate nodes, local index paired with the world + value (SC-007). Orthogonality rows with `f=y'`, direction `(1,0,0)`, + target 0, and direction `(0,1,0)`, target 1: residual at most `1e-10` + (SC-008, FDI and PLI; barycentre, `"centroid"` and `"gauss4"`). Feature + inequality bounds equal `h` at each node's world location. Restricted + reference support: `excluded_rows` count equals the analytic count of + samples outside the closed restricted box. Support expansion through a + restoring fault: restored world samples project inside the support + (FR-009, U2.7). + - `tests/unit/modelling/test_fault_frame_analytic.py`: R2 through + `GeologicalModel.create_and_add_fault` for FDI and PLI, explicit and + inferred centre, 4 anchors x 3 rotations. SC-003 on the 6859 probes with + CG `rtol=1e-10`, `atol=1e-12`, `maxiter=20000`. U1.1: generated rows use + the centre with no projection. Inferred centre equals `O` in all 3 + components. R4: the extent-direction row equals the analytic extent + gradient at the centre within `1e-10` per component; the analytic extent + directions at `O+(0,-6,0)` and `O+(0,6,0)` differ by at least 15 degrees; + a constant centre direction at those probes fails the check (SC-013). + Non-orthogonal slip input: `fault_slip_vector` is orthogonal to the + normal; a parallel slip raises `ValueError`. + - `tests/unit/modelling/test_fault_buffer_bounds.py`: R5 for FDI and PLI, + buffers 0, 0.2, 0.5, 3 rotations, explicit and inferred centre, forced + geometry. Each local face is outward of the expected enclosure by at most + 1 cell plus `Tp`. The 3 buffers give distinct bounds. All generated rows + are inside the support (SC-005). A property change after creation + (`fault_builder.fault_normal_vector = ...`) keeps the buffer and the + forced geometry (U1.5). + - `tests/fixtures/coordinate_frames.py` (new helpers, no tests): anchors, + `rot_z(deg) -> np.ndarray`, `tp(L, W) -> float`, R1 to R5 data frames, + the R2 analytic frame, the R3 field `h`, and the R4 curved fields. Tests + import only this module, not an external investigation workspace. +- **Integration and compatibility tests**: `tests/unit/test_public_api_contract.py` + and `tests/unit/test_stable_api_surface.py` must pass unchanged. If + `tests/fixtures/api_surface_snapshot.json` lists public methods of + `GeologicalInterpolator` or `BoundingBox`, add the new names to it (additive). +- **Round-trip tests**: S1 covers point and direction round trips (SC-001). + Persistence round trips are S3. +- **Backend coverage**: S1 covers Python FDI and PLI/P1 cells for R2, R3 and + R5. S2 and S6 cover the other cells (OP1). +- **Cross-language fixtures**: Not applicable in S1. +- **Grid and display checks**: S1 checks support enclosure and resize. Cell + centres and display are S4 and S5. +- **Performance check**: none required. Record the S1 test time (OP2). +- **Merge gates**: `pytest` and `ruff check LoopStructural packages`. Also run + `pytest packages/loop_common/tests packages/loop_interpolation/tests`, + because the root `testpaths` is `tests` only. + +## Project structure + +### Common geometry (`packages/loop_common`) + +**Responsibility**: rigid-rotation meaning, transformed box enclosure, stable +support resize. + +**Files**: + +- `packages/loop_common/src/loop_common/geometry/_rigid.py` - Create; + `RigidRotationError`, `validate_rigid_rotation`. +- `packages/loop_common/src/loop_common/geometry/__init__.py` - Change; export + `RigidRotationError` and `validate_rigid_rotation`. +- `packages/loop_common/src/loop_common/geometry/_bounding_box.py` - Change; + validation in `set_local_transform`; new `local_enclosure`. +- `packages/loop_common/src/loop_common/geometry/_structured_grid_3d.py` - + Change; `rotation_xy` validation; `_cells_for_length`; `origin` and + `maximum` setters. +- `packages/loop_common/tests/test_rigid_transform_contract.py` - Create; + SC-001, SC-011 (resize part), SC-012 (rotation part), enclosure. + +### Interpolation (`packages/loop_interpolation`) + +**Responsibility**: public world/local conversion; documented local and world +callback contracts. + +**Files**: + +- `packages/loop_interpolation/src/loop_interpolation/_geological_interpolator.py` + - Change; four public conversion methods. +- `packages/loop_interpolation/src/loop_interpolation/_discrete_interpolator.py` + - Change; `add_inequality_feature` samples world nodes; contract docstring + of `add_gradient_orthogonal_constraints`. + +### Geological modelling (`LoopStructural`) + +**Responsibility**: world fault geometry, reference sampling, support +requests. + +**Files**: + +- `LoopStructural/modelling/features/builders/_geological_feature_builder.py` + - Change; `excluded_rows`; `install_gradient_constraint`; + `install_equality_constraints`; `set_interpolation_geometry`; + `check_interpolation_geometry`. +- `LoopStructural/modelling/features/builders/_fault_builder.py` - Change; + stored buffer and forced geometry; 3D inferred centre; slip + orthogonalisation; extent-direction row; early-return mesh geometry. +- `LoopStructural/modelling/core/_model_feature_factory.py` - Change; remove + the centre projection and the buffer override in `create_and_add_fault`. +- `tests/fixtures/coordinate_frames.py` - Create; shared R1 to R5 fixtures. +- `tests/unit/modelling/test_reference_constraint_frames.py` - Create; R3, + FR-004 to FR-007, FR-009. +- `tests/unit/modelling/test_fault_frame_analytic.py` - Create; R2, R4, + FR-010, FR-011, SC-003, SC-013. +- `tests/unit/modelling/test_fault_buffer_bounds.py` - Create; R5, FR-012, + FR-013, SC-005. +- `tests/unit/modelling/test__fault_builder.py` - Change only if an assertion + encodes a defect fixed by S1 (for example an inferred centre with `z=0`). + Record each change in `verification.md` under C-004. + +### Documents + +**Responsibility**: evidence, golden review, migration notes. + +**Files**: + +- `specs/005-coordinate-frame-correctness/verification.md` - Create; + sections "Pre-fix failures", "Post-fix results", "Golden review (C-004)", + "Backend matrix", "Rollback checks", "Test run time (OP2)". +- `COMPAT.md` - Change; one entry under "Migration notices": explicit + `fault_center` is a world point; `fault_buffer` and `force_mesh_geometry` + now keep their effect; an inferred centre keeps its elevation; non-rigid + rotations are rejected; extra observations added as a workaround are no + longer necessary. + +### Later stages (file survey only, OP1) + +The next planner runs read and confirm these files. The lists come from a +grep survey on 2026-11-04 and are not final. + +- **S2**: `packages/loop_interpolation/src/loop_interpolation/_fd_fold_interpolator.py`, + `_discrete_fold_interpolator.py`, `_p2_fold_interpolator.py`, + `_p2interpolator.py`, `_fold_setup.py`, `_fold_norm_alignment.py`, + `_quadrature.py`, `_regularisation.py`, `_constant_norm.py`, + `_finite_difference_interpolator.py`, `_p1interpolator.py`; + `LoopStructural/modelling/features/builders/_folded_feature_builder.py`. +- **S3**: `_geological_interpolator.py` (`to_dict`, `from_dict`), + `_interpolator_factory.py`, `_finite_difference_interpolator.py` (`copy`), + `_p2interpolator.py` (`copy`), `_p2_fold_interpolator.py` (`copy`), + `packages/loop_common/src/loop_common/supports/_support_factory.py`, + `LoopStructural/modelling/features/fault/_fault_segment.py` (restoration + cache), `ROADMAP.md`, `COMPAT.md`. +- **S4**: `packages/loopstructural_visualisation/src/loopstructuralvisualisation/_3d_viewer.py`, + `LoopStructural/modelling/intrusions/intrusion_builder.py`, + `LoopStructural/modelling/intrusions/intrusion_frame_builder.py`, + `LoopStructural/modelling/features/fault/_fault_segment.py`, + `LoopStructural/export/exporters.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`. +- **S6**: `packages/loop_common/src/loop_common/math/_transformation.py`, + `LoopStructural/utils/_transformation.py`, + `packages/loop_interpolation/src/loop_interpolation/_native.py`, + `_surfe_wrapper.py`, `LoopStructural/modelling/features/fault/_fault_function.py`, + `ROADMAP.md`. + +## Requirement coverage + +| Requirement | Design and files | Verification | Gap | +|---|---|---|---| +| FR-001 | Existing public evaluation; S1 keeps it; fault path in `_model_feature_factory.py` | `test_fault_frame_analytic.py` (gradients in world after rotation) | None for S1 paths; other paths OP1 | +| FR-002 | Factory keeps explicit centre in world; builder rows in world | `test_fault_frame_analytic.py` U1.1 | None | +| FR-003 | Public conversion methods; builder uses each once | `test_rigid_transform_contract.py`; `test_reference_constraint_frames.py` | Load path (U3.2) is S3, OP1 | +| FR-004 | `install_gradient_constraint`, `install_equality_constraints`, `add_inequality_feature`; documented callbacks | `test_reference_constraint_frames.py` | None | +| FR-005 | `world_to_local_vectors` in `install_gradient_constraint` | SC-007, SC-008 tests | P2 is S2, OP1 | +| FR-006 | `install_equality_constraints` | 168-node check, rotated sets | None | +| FR-007 | `add_inequality_feature` | Bounds equal `h` at world nodes | None | +| FR-008 | S2 | S2 | OP1 | +| FR-009 | `check_interpolation_geometry` | U2.7 test | None | +| FR-010 | Factory (no `model.scale`); 3D inferred centre | SC-003 explicit and inferred | None | +| 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-016 | S3 | S3 | OP1 | +| FR-017 | S3 | S3 | OP1 | +| FR-018 | S3 | S3 | OP1 | +| FR-019 | S3 | S3 | OP1 | +| FR-020 | S3 | S3 | OP1 | +| FR-021 | `validate_rigid_rotation` in `set_local_transform` and `rotation_xy` | `test_rigid_transform_contract.py` | Saved-metadata part S3, OP1 | +| FR-022 | S4 | S4 | OP1 | +| FR-023 | S5 | S5 | OP1 | +| FR-024 | `excluded_rows` (S1 part) | Restricted-support exclusion count | Report integration S5, OP1 | +| 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 | +| FR-028 | No signature change; `COMPAT.md` note | API contract tests | FR-019 exception S3, OP1 | +| FR-029 | `verification.md` golden review section | Review record per changed result | None | +| FR-030 | S1 FDI and PLI cells | S1 tests | Other cells S2, S6, OP1 | +| FR-031 | Pre-fix rule in "Test strategy" | `verification.md` | None | +| FR-032 | "Stages and rollback" | Rollback checks in `verification.md` | None | +| SC-001 | `_rigid.py`, `BoundingBox` | `test_rigid_transform_contract.py` | None | +| SC-002 | R1 evaluation | S2/S3 (R1 lifecycle and fits) | OP1 | +| SC-003 | Fault design | `test_fault_frame_analytic.py` | Run time OP2 | +| SC-004 | S6 | S6 | OP1 | +| SC-005 | Fault and enclosure design | `test_fault_buffer_bounds.py` | Run time OP2 | +| SC-006 | S3 | S3 | OP1 | +| SC-007 | Reference sampling | `test_reference_constraint_frames.py` | Fold and P2 S2, OP1 | +| SC-008 | Reference sampling | `test_reference_constraint_frames.py` (FDI, PLI) | P2 S2, OP1 | +| SC-009 | S5 | S5 | OP1 | +| SC-010 | S6 | S6 | OP1 | +| SC-011 | `_cells_for_length` (resize part) | Resize test | Cell centres S5, OP1 | +| SC-012 | Rotation part in S1 | `test_rigid_transform_contract.py` | Metadata and legacy part S3, OP1 | +| SC-013 | Extent-direction row | `test_fault_frame_analytic.py` | Fold R4 geometry S2, OP1 | +| SC-014 | Merge gates in every stage | `pytest`, ruff, package tests | None | +| SC-015 | `verification.md` golden review | Review record; 0 tolerance changes | None | +| SC-016 | S3 | S3 | OP1 | +| SC-017 | Rollback column; format part S3 | `verification.md` | OP1 | +| SC-018 | S1 cells; matrix S6 | `verification.md` backend matrix | OP1 | +| SC-019 | Gate checked on 2026-11-04: 003 spec (C-006, main and native worktree copies identical) and 004 spec (2026-11-03 source-005 entries) contain the amendments | Document check | None | + +## Handover + +- **Coordinate conversion (owner: `loop_interpolation`)**: + `GeologicalInterpolator.world_to_local_points`, `local_to_world_points`, + `world_to_local_vectors`, `local_to_world_vectors`. Spec 004 uses + `world_to_local_vectors` to convert sampled world directions once into the + computational basis before the low-level core. Spec 001 keeps its contract: + the core gets local, unrotated grid inputs. +- **Rotation check (owner: `loop_common`)**: + `loop_common.geometry.validate_rigid_rotation(matrix, dimensions, atol=1e-10)` + and `RigidRotationError(ValueError)`. S3 uses them for saved metadata. +- **Box enclosure (owner: `loop_common`)**: + `BoundingBox.local_enclosure(origin=None, maximum=None) -> (local_min, local_max)`. +- **Callback contracts**: feature-facing reference callbacks and region + callbacks receive world points. `add_gradient_orthogonal_constraints` + receives target-local points and target-local unit vectors. +- **Exclusion record**: `GeologicalFeatureBuilder.excluded_rows: + dict[str, dict[str, int]]`. S5 connects it to the diagnostics report. +- **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 + record in `verification.md`. Spec 003 cache invalidation for transform + changes is S3. +- **Shared files**: `_geological_interpolator.py`, `_discrete_interpolator.py` + and `_fault_builder.py` are also in the scope of specs 003 and 004. + Coordinate the merge order with the active 003 worktree + (`LoopStructural-native`). + +## Complexity tracking + +None. diff --git a/specs/005-coordinate-frame-correctness/run-state.json b/specs/005-coordinate-frame-correctness/run-state.json new file mode 100644 index 000000000..856b2b13a --- /dev/null +++ b/specs/005-coordinate-frame-correctness/run-state.json @@ -0,0 +1,351 @@ +{ + "spec": "005-coordinate-frame-correctness", + "stage": "S1", + "worktree": "C:/Users/lachl/dev/LoopStructural-005", + "branch": "feat/005-coordinate-frame-correctness", + "review_mode": "wave", + "start_commit": "3881ef320adc483f71326ef9d09d90596493d39f", + "baseline": { + "lint": "pass: ruff check LoopStructural packages", + "pytest": "755 passed, 5 skipped, 21 warnings in 1205.98 s", + "doctor": "pass: source CLI workaround from setup-notes.md" + }, + "runtime": { + "max_workers": 3, + "worker_S": "gpt-6-luna low", + "worker_M": "gpt-6-luna medium", + "review_wave": "gpt-6-luna medium", + "review_spec": "gpt-6.1-sol medium" + }, + "waves": [ + { + "wave": 1, + "base": "3881ef320adc483f71326ef9d09d90596493d39f", + "tasks": [ + "T001", + "T002", + "T003" + ], + "review": "pass: gpt-6-luna medium, no findings", + "check": "scope pass; verification skeleton pass; fixture 0 tests; transform import fails as required before T004; lint and format pass", + "commit": "51c7c1a3553a302dfc0add30a1b97434e9645761" + }, + { + "wave": 2, + "base": "51c7c1a3553a302dfc0add30a1b97434e9645761", + "tasks": [ + "T004", + "T007", + "T011", + "T014" + ], + "review": "pass after one fixed major finding; fix recheck gpt-6-luna medium: no findings", + "check": "scope, Black and Ruff pass; validator 13 pass; parent regression run 373 expected failures, 28 pass, 48 analytic cells deselected, 39.85 s; worker previously ran all 48 analytic cells red", + "commit": "dd7078d62b7eedf8cde4e453ac622fda8ccafb73", + "fix_commit": "72d5e1b8dec1e78309cafc9258b1e761bad03783" + }, + { + "wave": 3, + "base": "72d5e1b8dec1e78309cafc9258b1e761bad03783", + "tasks": [ + "T005", + "T006", + "T008", + "T012" + ], + "review": "pass: sonnet, 1 minor (T012 asarray in-place change of caller normal)", + "status": "committed", + "check": "scope pass (T006 golden file added to Files by user decision); ruff pass; black pass; loop_common + loop_interpolation tests 912 passed, 47 skipped, 2 xfailed; conversion filter 24 pass; T012 filter 4 pass; test_base.py errors not reproduced", + "commit": "d2909aa5de619d9e83c49b818dfc7d1333700d51" + }, + { + "wave": 4, + "base": "d2909aa5de619d9e83c49b818dfc7d1333700d51", + "tasks": [ + "T009", + "T010" + ], + "check": "scope pass; ruff pass; tests/unit/modelling (without T013 files) 551 passed, 2 skipped; reference file 312 pass; loop_interpolation tests 606 passed, 47 skipped, 2 xfailed", + "review": "pass: sonnet, no findings (2 notes, not findings)", + "commit": "35e05b01c858f619cb667af446f507046cc78fb3" + }, + { + "wave": 5, + "base": "35e05b01c858f619cb667af446f507046cc78fb3", + "tasks": [ + "T013" + ], + "check": "scope pass; ruff pass; tests/unit/modelling 688 passed, 2 skipped in 300 s; R2+R5 files 137 passed in 209 s", + "review": "pass: sonnet, no findings", + "commit": "6d533fe01315647780a7e22c7c6bed2e09d9c076" + }, + { + "wave": 6, + "base": "6d533fe01315647780a7e22c7c6bed2e09d9c076", + "tasks": [ + "T015", + "T016" + ], + "check": "scope pass; test__fault_builder.py 6 pass (unchanged, no golden change); golden npz entry matches orchestrator key comparison; COMPAT row in Migration notices table", + "review": "pass: sonnet, 1 minor (golden entry omits 33 round-off keys) - fixed in T017", + "commit": "10fd4dea59d47c86da99aad90c66e437092ffe3b" + }, + { + "wave": 7, + "base": "10fd4dea59d47c86da99aad90c66e437092ffe3b", + "tasks": [ + "T017" + ], + "check": "root ruff pass; root pytest 1204 passed, 5 skipped in 1082.66 s (orchestrator, HEAD 10fd4dea); packages 912 passed, 47 skipped, 2 xfailed; S1 files 491 passed in 291 s; rollback at ebc7c699 15 passed; all 6 verification headings filled; prefix worktree removed", + "commit": "0abc31ecc4d95d1a2b7c205701eb265416c23a2c" + } + ], + "runs": { + "T001": [ + { + "model": "gpt-6-luna", + "reasoning_effort": "low", + "result": "pass: verification skeleton read back" + } + ], + "T002": [ + { + "model": "gpt-6-luna", + "reasoning_effort": "medium", + "result": "pass: imports, 0 tests collected, lint and format pass" + } + ], + "T003": [ + { + "model": "gpt-6-luna", + "reasoning_effort": "medium", + "result": "needs changes: missing 2D/rotation cases and exact vector tolerance" + }, + { + "model": "gpt-6-luna", + "reasoning_effort": "medium", + "result": "pass: expected missing RigidRotationError import; lint and format pass" + }, + { + "model": "sonnet", + "result": "fix pass: dimensions=dimensions in round-trip BoundingBox; contract file 42 pass" + } + ], + "T004": [ + { + "model": "gpt-6-luna", + "reasoning_effort": "low", + "result": "needs changes: targeted 13 pass; Google docstring and lint suppression must be corrected" + }, + { + "model": "gpt-6-luna", + "reasoning_effort": "low", + "result": "pass: NumPy docstring; sorted imports without suppression; orchestrator verified 13 tests, lint pass" + } + ], + "T007": [ + { + "model": "gpt-6-luna", + "reasoning_effort": "medium", + "result": "fail: guessed constraint internals, incomplete matrix and weak assertions" + }, + { + "model": "gpt-6-luna", + "reasoning_effort": "medium", + "result": "pass: corrected real APIs, complete Cartesian matrix, fixed world affine controls, wrapper checks, physical restoration, and exclusion counts; lint/format pass" + }, + { + "model": "sonnet", + "result": "fix pass: exact field without large-coordinate cancellation, atol 1e-10 kept; file 312 pass" + } + ], + "T014": [ + { + "model": "gpt-6-luna", + "reasoning_effort": "medium", + "result": "needs changes: 42 expected pre-fix failures but PLI and generated-row containment coverage missing" + }, + { + "model": "gpt-6-luna", + "reasoning_effort": "medium", + "result": "pass: both engines, all supports, world-row containment, correct Tp, exact unchanged state and py39; parent independently confirmed intended failures" + } + ], + "T011": [ + { + "model": "gpt-6-luna", + "reasoning_effort": "medium", + "result": "needs changes: 52 expected failures, 1 pass, all 48 analytic cells ran in 145.29 s; missing nonorthogonal input and finite-probe guards" + }, + { + "model": "gpt-6-luna", + "reasoning_effort": "medium", + "result": "pass: actual nonorthogonal input, single extent row, all-finite probes; 4 expected representative failures; Black/Ruff pass" + }, + { + "model": "gpt-6-luna", + "reasoning_effort": "medium", + "result": "review fix pass: extent-row test uses inferred centre to isolate T012; parent confirms expected zero-row failure, lint pass" + }, + { + "model": "sonnet", + "result": "fix pass: slip row expectation [0,0,-1]/intermediate_axis per spec.md R2 table; T012 filter 4 pass" + }, + { + "model": "sonnet", + "result": "fix pass: SC-003 CG settings through solver_kwargs (top-level rtol/atol/maxiter were ignored; solver used atol 2e-9); R2+R5 137 pass" + } + ], + "T005": [ + { + "model": "gpt-6-luna", + "reasoning_effort": "medium", + "result": "stopped, edits saved: enclosure/prior-state 4 pass; 12 round-trip failures from missing dimensions=2 in T003 test; Black/Ruff pass; common package incomplete with 3 reported test_base.py errors, no traceback retained" + } + ], + "T006": [ + { + "model": "gpt-6-luna", + "reasoning_effort": "medium", + "result": "stopped, edits saved: prior_state/unchanged_bounds 2 pass, 40 deselected; Black pass; full common package cancelled with Ctrl+C at user request, exit 1, no final result" + }, + { + "model": "sonnet", + "result": "fix pass: ruff I001; golden npz regenerated, 27 keys changed (after_origin/after_maximum only); common tests 306 pass" + } + ], + "T012": [ + { + "model": "gpt-6-luna", + "reasoning_effort": "medium", + "result": "stopped, edits saved: focused inferred_centre/slip/extent_direction 4 pass, 49 deselected before formatting; Black/Ruff pass; parent verification pending" + }, + { + "model": "sonnet", + "result": "fix: restored coord-1 slip gradient scale 1/intermediate_axis (spec R2 table); blocked on wrong T011 expectation" + } + ], + "T008": [ + { + "model": "haiku", + "result": "pass: 4 public conversion methods; conversion filter 24 pass" + } + ], + "T009": [ + { + "model": "haiku", + "result": "pass: inequality 48 pass; task hunks only, no full-file black reformat" + } + ], + "T010": [ + { + "model": "sonnet", + "result": "blocked: 12 affine-residual failures at anchor 6.5e6 from test cancellation; code otherwise complete" + } + ], + "T013": [ + { + "model": "haiku", + "result": "code as Context; Done-when failed 48 R2 cells at 1e-4 because the T011 test passed CG settings outside solver_kwargs; pass after T011 fix" + } + ], + "T015": [ + { + "model": "sonnet", + "result": "pass: no assertion in test__fault_builder.py encodes an S1 defect; npz golden entry recorded" + } + ], + "T016": [ + { + "model": "haiku", + "result": "pass: row added at end of Migration notices table" + } + ], + "T017": [ + { + "model": "sonnet", + "result": "pass: pre-fix evidence (contract not runnable, 278/53/84 failures), rollback 15 pass, timing, backend matrix, golden round-off note" + } + ] + }, + "findings": { + "found": 13, + "rejected": 0, + "dropped": 0, + "fixed": 6, + "open_minor": [ + "T011 test_extent_direction_row_at_centre: 15-degree check uses only the r4_frame fixture; production row not compared at O+(0,+-6,0) (SC-013 negative case)", + "T011 test_explicit_centre_rows_not_projected: builders[0] check is always true (R2 has an observation at the centre)", + "T007 test_equality_region_selects_world_nodes: 168 of 336 count for rotation 0 not asserted", + "T007 test_set_interpolation_geometry_encloses_rotated_box: passes rotation in degrees, which also sets support.rotation_xy (double rotation in setup)", + "T010/T005: no test for the evaluation_error branch of excluded_rows or for non-finite local_origin rejection (code checked correct by scratch test)", + "Later stage note: _on_geometry_change passes fault_trace_anisotropy=1.0 and fault_dip_anisotropy=1.0, so a property change adds anisotropy rows even when the caller set 0", + "T011 test_caller_normal_not_mutated: the slip assertion cannot fail, because the factory copies slip first; only the normal assertion checks the fix" + ] + }, + "implementation_notes": [ + "T004 import position follows project Ruff ordering without suppression; validator behavior and exports match the plan.", + "Orchestrator diagnosis wave 5: with solver_kwargs the anchor-0 R2 error is 2e-12 to 1e-10 (FDI and PLI); ridge regularisation has no effect. Top-level rtol/atol/maxiter on builder.build are ignored silently; possible later issue, not in S1 scope." + ], + "pre_fix_evidence": { + "T011": [ + "inferred centre [100,200,0] versus [100,200,300]", + "explicit generated row [3.66025404,13.66025404,10] versus world [100,200,300]", + "nonorthogonal slip [0.7071,0,-0.7071] instead of [0,0,-1]", + "zero finite coord2 direction rows at centre", + "all 48 analytic matrix cells fail before fixes" + ], + "T014": [ + "b=0 actual support origin [-4,-4,-4] versus requested [0,0,0], spacing [1,1,1], Tp 4.82e-8", + "90-degree frames raise resize ValueError before fixes" + ], + "wave2_parent": "373 expected failures, 28 pass, 48 deselected in 39.85 s; first failure of each reference, centre, slip, extent and buffer group independently confirmed" + }, + "review_fixes": [ + { + "wave": 2, + "task": "T011", + "severity": "major", + "problem": "extent test depended on later explicit-centre factory fix", + "fix": "use inferred centre for isolated extent test", + "status": "fixed and recheck pass" + } + ], + "runtime_resume": { + "date": "2026-11-04", + "runtime": "claude", + "max_workers": 5, + "worker_S": "haiku", + "worker_M": "sonnet", + "review_wave": "sonnet", + "review_spec": "opus", + "note": "Resumed in Claude Code after the Codex handoff; models follow providers.claude of .agents/model-policy.json." + }, + "answers": [ + "2026-11-04: golden structured_support_reference.npz added to T006 Files; T015 records it (plan.md Answers)" + ], + "spec_review": { + "model": "opus", + "range": "3881ef32..0abc31ec", + "verdict": "changes needed", + "major": 2, + "minor": 7, + "fixed": [ + "major: T012 caller normal changed in place (np.array copy) + test_caller_normal_not_mutated", + "major: T013 fault_ellipsoid double conversion of world centre (_fault_segment.py added to T013 Files by user) + test_fault_ellipsoid_centre_is_world", + "minor (user asked): verification.md Pre-fix titles FR ranges corrected", + "minor (user asked): golden npz keeps only the 27 FR-014 changes; 33 round-off keys restored" + ], + "fix_check": "scope pass; ruff pass; tests/unit/modelling + loop_common (without planar matrix) 948 passed, 2 skipped", + "commit": "6547516e0090b311598fd94d7df73e2fe590881a", + "recheck": "pass: opus, 1 minor (slip assertion in test_caller_normal_not_mutated cannot fail; factory already copies slip)" + }, + "final": { + "head": "6547516e0090b311598fd94d7df73e2fe590881a", + "lint": "ruff check LoopStructural packages: pass", + "format": "black --check on the 15 changed .py files: 13 unchanged, 2 would be reformatted (_fault_segment.py, _discrete_interpolator.py; both not black-clean before S1, workers kept task hunks only)", + "test_root": "pytest: 1206 passed, 5 skipped, 21 warnings in 867 s (baseline 755 passed, 5 skipped in 1206 s)", + "test_packages": "pytest packages/loop_common/tests packages/loop_interpolation/tests: 912 passed, 47 skipped, 2 xfailed in 117 s", + "status": "S1 implemented; S2-S6 need planning (OP1)" + } +} diff --git a/specs/005-coordinate-frame-correctness/spec.md b/specs/005-coordinate-frame-correctness/spec.md index 770e41c00..40bbe4adf 100644 --- a/specs/005-coordinate-frame-correctness/spec.md +++ b/specs/005-coordinate-frame-correctness/spec.md @@ -1,6 +1,6 @@ # Feature spec: Coordinate and fault-frame correctness -**Feature**: `005-coordinate-frame-correctness` | **Created**: 2026-11-03 | **Status**: Review corrections applied; planning blocked on the required cascade to specs 003 and 004 +**Feature**: `005-coordinate-frame-correctness` | **Created**: 2026-11-03 | **Status**: Stage S1 implemented (pull request #392); stages S2 to S6 need planning **Input**: Fix the shared world-to-local coordinate defects and the fault-frame interpolation defects found in the two investigations. @@ -402,6 +402,8 @@ An unsupported capability must retain its documented rejection. Missing optional - A new RBF, native, browser, or GPU solver; performance work beyond checking that correctness changes do not remove existing supported paths. - A rewrite of the intrusion workflow, QGIS plugin, model graph, or all export formats. The audited coordinate consumers and shared geometry are in scope. - Changes to other spec documents, branches, worktrees, commits, or release publication during this interview. +- The trace and dip anisotropy settings of a fault after a property change. A rebuild after a property change uses anisotropy weights of 1.0, also when the caller set 0. FR-012 covers only the fault buffer and forced mesh geometry. Issue #394 records this defect (C-005). +- Solver settings that `build` ignores. Top-level `rtol`, `atol` and `maxiter` have no effect, and no warning occurs; only `solver_kwargs` applies them. This is not a coordinate-frame defect. Issue #393 records it (C-006). The tests of this spec give the SC-003 solve settings through `solver_kwargs`. ## Risks and rollback @@ -439,6 +441,14 @@ All four interview decisions are resolved. The requirements, acceptance scenario - Made the accepted cascade to specs 003 and 004 a pre-plan prerequisite. This review has not changed those documents. - Review verdict remains **Needs changes** only because that external prerequisite is still outstanding. No new user-policy clarification is required. +### Decisions after stage S1, answered on 2026-11-04 + +The S1 build (pull request #392) found two defects that the spec did not decide, and the maintainer decided the plan point OP2. + +- **C-005, anisotropy after a property change**. Question: must FR-012 also keep the trace and dip anisotropy settings through a rebuild? Answer: no. This is out of scope for spec 005. Issue #394 records it. +- **C-006, ignored solver settings**. Question: must this spec make `build` apply or reject top-level solver settings? Answer: no. It is not a coordinate-frame defect. Issue #393 records it. +- **C-007, R2 and R5 run time (plan OP2)**. Question: where does the mandatory R2 and R5 matrix run? Answer: in the default `pytest` run, with no marker. In S1, the four new test files take about 5 minutes. The analytic file takes 264 s. The full root suite takes about 14.5 minutes. + ## Next step -Run spec-cascade to apply the accepted changes to specs 003 and 004. Recheck their consistency with 005 and confirm SC-019. Planning may start only after that gate passes. No clarification marks remain. +Stage S1 is implemented in pull request #392. After it merges, run the spec-planner skill for stage S2 (plan OP1). No clarification marks remain. diff --git a/specs/005-coordinate-frame-correctness/tasks.md b/specs/005-coordinate-frame-correctness/tasks.md new file mode 100644 index 000000000..ad9e750ff --- /dev/null +++ b/specs/005-coordinate-frame-correctness/tasks.md @@ -0,0 +1,403 @@ +# Tasks: Coordinate and fault-frame correctness (stage S1) + +**Feature**: `005-coordinate-frame-correctness` | **Plan**: [plan.md](plan.md) + + + +## 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: 7. Tiers: [S] = simple model, [M] = standard model. +- Do not start a task with "Blocked by". + +## Phase 1: Setup + +- [x] T001 [P] [Setup] [S] Create the S1 verification record skeleton. + - Files: `specs/005-coordinate-frame-correctness/verification.md` + - Needs: none + - Wave: 1 + - Done when: the file exists with the 6 headings below, each followed by + the line "No entries yet." + - Context: Title `# Verification: 005-coordinate-frame-correctness`. + Headings, in this order: `## Pre-fix failures`, `## Post-fix results`, + `## Golden review (C-004)`, `## Backend matrix`, `## Rollback checks`, + `## Test run time (OP2)`. Under the title, add one sentence: "This file + records the evidence for FR-029, FR-031 and FR-032 of spec 005." + Write in ASD-STE100 Simplified Technical English. + +- [x] T002 [P] [Setup] [M] Create the shared analytic fixtures R1 to R5. + - Files: `tests/fixtures/coordinate_frames.py` + - Needs: none + - Wave: 1 + - Done when: `pytest tests/fixtures/coordinate_frames.py` collects 0 + tests with no import error, and `ruff check LoopStructural packages` + passes. + - Context: Helper module, no tests. numpy docstrings and type hints; Black + line length 100. Names: `ANCHORS` = arrays (0,0,0), (100,200,300), + (-100,-200,-300), (500000,6500000,100); `ROTATIONS_DEG = (0.0, 30.0, 90.0)`; + `rot_z(deg: float) -> np.ndarray` = [[c,-s,0],[s,c,0],[0,0,1]]; + `tp(L: float, W: float) -> float` = `1e-10*L + 32*np.finfo(float).eps*max(W, L)`; + `model_with_frame(origin, maximum, rotation_deg) -> GeologicalModel` makes + `GeologicalModel(origin, maximum)` then calls + `model.bounding_box.set_local_transform(local_origin=origin, rotation_matrix=rot_z(rotation_deg))`. + `r2_fault_data(anchor)`: DataFrame, columns X,Y,Z,val,coord,feature_name + ("fault"),w; 25 rows at x'=0, y',z' in {-8,-4,0,4,8}, val 0, coord 0, w 1. + `r5_fault_data(anchor)`: same, y',z' in {-4,-2,0,2,4}. + `r2_probes(anchor)`: 6859 points `anchor + {-9..9}^3`. + `r2_frame(points, anchor) -> (values (N,3), gradients (N,3,3))`: + gx=x'/20, gy=-z'/20, gz=y'/20. `r3_value`/`r3_gradient`: h=x'+0.05*y'^2+z', + grad (1,0.1*y',1). `r4_frame(points, anchor)`: R=40, X=R+x', + r=sqrt(X^2+y'^2); surface (r-R)/20, slip -z'/20, extent R*atan2(y',X)/20, + with the analytic gradients of spec.md R4. `class AnalyticReference` with + `__init__(value_fn, gradient_fn, name, support_box=None)`, + `evaluate_value(points)`, `evaluate_gradient(points, ignore_regions=True)`, + `__call__ = evaluate_value`; NaN outside the closed `support_box`. + +## Phase 2: Foundation (blocks all user stories) + +- [x] T003 [P] [Found] [M] Test: rigid rotation, round trips, enclosure and resize stability. + - Files: `packages/loop_common/tests/test_rigid_transform_contract.py` + - Needs: none + - Wave: 1 + - Done when: `pytest packages/loop_common/tests/test_rigid_transform_contract.py` + fails before T004 (import error), then passes after T005 and T006. + - Context: Self-contained (do not import `tests/fixtures`). Tests: + `test_validate_rigid_rotation_accepts_proper_rotations` (0, 30, 90 deg + about z; 2D and 3D); `test_validate_rigid_rotation_rejects_invalid` + (scale, shear, singular, reflection det -1, NaN, inf, wrong shape) expects + `RigidRotationError`, a `ValueError` subclass; + `test_point_and_direction_round_trips` (SC-001: `BoundingBox.project` then + `reproject` within `Tp = 1e-10*L + 32*eps*max(W,L)`; `project_vectors` then + `reproject_vectors` within 1e-12; 1 point and arrays; anchors (0,0,0), + (100,200,300), (-100,-200,-300), (500000,6500000,100)); + `test_local_enclosure_contains_projected_corners` (all 8 projected corners + inside `BoundingBox.local_enclosure()` within Tp, for each rotation); + `test_invalid_rotation_keeps_prior_state` (`set_local_transform` and + `StructuredGrid3DGeometry.rotation_xy` keep the old matrix after the error); + `test_reapply_unchanged_bounds_is_stable` (set `origin` and `maximum` of a + `StructuredGrid3DGeometry` to their own values 10 times: 0 changes to + origin, maximum, step_vector, nsteps). Import from `loop_common.geometry`. + +- [x] T004 [Found] [S] Add `validate_rigid_rotation` and `RigidRotationError`. + - Files: `packages/loop_common/src/loop_common/geometry/_rigid.py`, + `packages/loop_common/src/loop_common/geometry/__init__.py` + - Needs: T003 + - Wave: 2 + - Done when: `pytest packages/loop_common/tests/test_rigid_transform_contract.py -k validate_rigid_rotation` + passes. + - Context: New module, imports numpy only. `class RigidRotationError(ValueError)`. + `def validate_rigid_rotation(matrix: np.ndarray, dimensions: int, atol: float = 1e-10) -> np.ndarray`: + return `np.array(matrix, dtype=float)` (a copy). Raise + `RigidRotationError` when shape != (dimensions, dimensions), any entry is + not finite, `max|R.T @ R - I| > atol`, or `|det(R) - 1| > atol`. Message + names the failed condition and ends with "supply a proper rotation + (orthonormal, det +1); scale and shear are not supported". In + `__init__.py` add `from ._rigid import RigidRotationError, validate_rigid_rotation` + as the first import line. numpy docstrings, type hints, py39 (use + `from __future__ import annotations`). + +- [x] T005 [P] [Found] [M] Validate the local transform and add `BoundingBox.local_enclosure`. + - Files: `packages/loop_common/src/loop_common/geometry/_bounding_box.py` + - Needs: T004 + - Wave: 3 + - Done when: `pytest packages/loop_common/tests/test_rigid_transform_contract.py -k "round_trips or local_enclosure or prior_state"` + and `pytest packages/loop_common/tests` pass. + - Context: In `set_local_transform(local_origin=None, rotation_matrix=None)` + call `validate_rigid_rotation(rotation_matrix, self.dimensions)` (import + from `._rigid`) and reject a non-finite `local_origin` with + `LoopValueError`, both before any field changes. New method + `def local_enclosure(self, origin: np.ndarray | None = None, maximum: np.ndarray | None = None) -> tuple[np.ndarray, np.ndarray]`: + defaults `self.origin`, `self.maximum`; raise `LoopValueError` if a bound + is None, not finite, or origin > maximum on an axis; build the + `2**dimensions` corners with `itertools.product`; return + `(min, max)` over axis 0 of `self.project(corners)`. numpy docstring. + Existing tests in `packages/loop_common/tests/test_bounding_box.py` use + proper rotations and must still pass. + +- [x] T006 [P] [Found] [M] Repair support resize and validate `rotation_xy`. + - Files: `packages/loop_common/src/loop_common/geometry/_structured_grid_3d.py`, + `packages/loop_common/tests/data/structured_support_reference.npz` + - Needs: T004 + - Wave: 3 + - Done when: `pytest packages/loop_common/tests/test_rigid_transform_contract.py -k "prior_state or unchanged_bounds"` + and `pytest packages/loop_common/tests` pass. + - Context: `_nsteps` is a NODE count; `nsteps_cells = nsteps - 1`; + `maximum = origin + nsteps_cells * step_vector`. Defect: the `origin` and + `maximum` setters set `_nsteps = ceil(length/step)` (a cell count), so the + support loses 1 cell. Add + `def _cells_for_length(self, length: np.ndarray) -> np.ndarray`: + `cells = np.ceil(length / self.step_vector - 1e-9).astype(np.int64)`; + axes with 0 cells get 2 cells; raise `ValueError("Cannot resize the grid.")` + if a count < 0. `origin` setter: compute `self.maximum` first, then + `self._nsteps = self._cells_for_length(maximum - origin) + 1`. + `maximum` setter: `self._nsteps = self._cells_for_length(maximum - self.origin) + 1`. + Keep `_on_geometry_change()` calls. `rotation_xy` setter: after the + float/int conversion call `validate_rigid_rotation(rotation_xy, 3)` + (import from `._rigid`); `None` keeps the current value. + Answer (plan.md "Answers", 2026-11-04): regenerate the golden file with + `python packages/loop_common/tests/test_structured_support_reference.py --update`; + only `after_origin`/`after_maximum` `nsteps` and `maximum` may change. + +## Phase 3: User story 2 - Constrain features at the intended physical points (P1) + +**Goal**: reference relationships sample at world points and install rows in +the target-local basis. +**Independent test**: `pytest tests/unit/modelling/test_reference_constraint_frames.py`. + +- [x] T007 [P] [US2] [M] Test: R3 reference relationships for FDI and PLI. + - Files: `tests/unit/modelling/test_reference_constraint_frames.py` + - Needs: T002 + - Wave: 2 + - Done when: `pytest tests/unit/modelling/test_reference_constraint_frames.py` + fails before T008 to T010, then passes after them. Record the failing + test ids in the task report. + - Context: Use `tests/fixtures/coordinate_frames.py` (ANCHORS, + ROTATIONS_DEG, model_with_frame, AnalyticReference, r3_value, r3_gradient, + tp). Box `O` to `O+(8,10,12)`; support 5x6x7 cells (336 nodes); engines + "FDI" and "PLI". Tests: `test_world_local_conversion_methods` + (`interpolator.world_to_local_points` etc. invert each other within Tp); + `test_equality_region_selects_world_nodes` (region `x' > 4`: 168 of 336 + unrotated; for rotations compute the expected set from world nodes; 0 + missing or duplicate; value = h at the world node); + `test_orthogonal_rows_affine_residual` (field f=y', direction (1,0,0) + target 0, direction (0,1,0) target 1; rows `f"{name}_orthogonal"`; + residual `A @ c_exact - b` at most 1e-10 for sample_rule None, + "centroid", "gauss4"); + `test_inequality_feature_bounds_world` (bounds = h at world nodes); + `test_restricted_reference_exclusion_counts` (box O+(2.1,2.1,2.1) to + O+(5.9,7.9,9.9); `builder.excluded_rows` count equals the analytic count); + `test_support_expansion_contains_restored_points` (FR-009); + `test_set_interpolation_geometry_encloses_rotated_box`. + +- [x] T008 [US2] [S] Add the public world/local conversion methods. + - Files: `packages/loop_interpolation/src/loop_interpolation/_geological_interpolator.py` + - Needs: T007 + - Wave: 3 + - Done when: `pytest tests/unit/modelling/test_reference_constraint_frames.py -k conversion` + passes. + - Context: Add four public methods after `_project_constraint_array`: + `world_to_local_points(self, points: np.ndarray) -> np.ndarray` returns + `self._project_points(points)`; `local_to_world_points` returns + `self._reproject_points(points)`; `world_to_local_vectors` returns + `self._project_vectors(vectors)`; `local_to_world_vectors` returns + `self._reproject_vectors(vectors)`. Keep the private methods. Docstrings + (numpy style): input `(N, d)` or `(d,)`, same output shape; vectors use + rotation only; with `bounding_box is None` the input is returned. + +- [x] T009 [P] [US2] [S] Sample `add_inequality_feature` at world nodes. + - Files: `packages/loop_interpolation/src/loop_interpolation/_discrete_interpolator.py` + - Needs: T008 + - Wave: 4 + - Done when: `pytest tests/unit/modelling/test_reference_constraint_frames.py -k inequality` + passes. + - Context: In `add_inequality_feature(self, feature, lower=True, mask=None)` + replace `value = feature(self.support.nodes)` with + `value = feature(self.local_to_world_points(self.support.nodes))`. Add to + its docstring: "`feature` receives world points and returns one value for + each support node." In the docstring of the abstract + `add_gradient_orthogonal_constraints(self, points, vectors, w=1.0)` state: + "`points` are target-local support points; `vectors` are target-local + unit directions." No signature changes. + +- [x] T010 [P] [US2] [M] Sample reference features at world points in the feature builder. + - Files: `LoopStructural/modelling/features/builders/_geological_feature_builder.py` + - Needs: T005, T008 + - Wave: 4 + - Done when: `pytest tests/unit/modelling/test_reference_constraint_frames.py -k "not inequality"` + and `pytest tests/unit/modelling` pass. + - Context: Add `self.excluded_rows: dict[str, dict[str, int]] = {}` in + `__init__`. Reasons: "invalid_reference_value", + "invalid_reference_direction", "evaluation_error". + `install_gradient_constraint`: local points (existing) -> `world = + self.interpolator.local_to_world_points(pts)`; `vector = + feature.evaluate_gradient(world, ignore_regions=True)`; `valid` = finite + and norm > 0; record `(~valid).sum()`; normalise; `world_to_local_vectors`; + pass `pts[valid]`, local vectors, `weight[valid]` if array. 0 valid -> no + rows, warning. `install_equality_constraints`: world nodes, `idc = + np.flatnonzero(region(world_nodes))`, `val = + feature.evaluate_value(world_nodes[idc])`, finite mask, name + `f"{feature.name}_equality"`; the `except` records "evaluation_error". + `set_interpolation_geometry`: use + `self.interpolator.bounding_box.local_enclosure(origin, maximum)`. + `check_interpolation_geometry`: `regular_grid(local=False)`, restore in + world, `world_to_local_points`, min/max on copies of support arrays. + +## Phase 4: User story 1 - Build a fault frame at its supplied location (P1) + +**Goal**: the public fault API makes the analytic planar frame at the supplied +world location, with the requested buffer. +**Independent test**: `pytest tests/unit/modelling/test_fault_frame_analytic.py tests/unit/modelling/test_fault_buffer_bounds.py`. + +- [x] T011 [P] [US1] [M] Test: R2 planar frame, centre handling and R4 extent row. + - Files: `tests/unit/modelling/test_fault_frame_analytic.py` + - Needs: T002 + - Wave: 2 + - Done when: `pytest tests/unit/modelling/test_fault_frame_analytic.py` + fails before T012 and T013, then passes after them. + - Context: Build through `model.create_and_add_fault("fault", 0, data=..., + interpolatortype=..., fault_center=O or None, fault_normal_vector=(1,0,0), + fault_slip_vector=(0,0,-1), major_axis=20, minor_axis=10, + intermediate_axis=20, fault_buffer=0.2, nelements=20000 (FDI) or 1000 + (PLI), fault_trace_anisotropy=0, fault_dip_anisotropy=0)`. Solve: CG rtol + 1e-10, atol 1e-12, maxiter 20000; orthogonality weight 1, step 10. + Tests: `test_explicit_centre_rows_not_projected` (generated rows use O); + `test_inferred_centre_keeps_elevation` (inferred centre == O, all 3 + components); `test_slip_vector_made_orthogonal` and + `test_parallel_slip_raises` (ValueError); + `test_extent_direction_row_at_centre` (one coord-2 row at O with + n = (0,1/20,0); equals the R4 analytic extent gradient at O within 1e-10; + R4 extent directions at O+(0,-6,0), O+(0,6,0) differ >= 15 degrees); + `test_planar_fault_frame_matches_analytic` (SC-003 on `r2_probes`: scalar + error <= 1e-6, relative gradient magnitude <= 1e-4, direction <= 0.01 deg; + 4 anchors x 3 rotations x 2 centre modes x FDI, PLI). + +- [x] T012 [US1] [M] Repair centre, slip, extent row and buffer state in `FaultBuilder`. + - Files: `LoopStructural/modelling/features/builders/_fault_builder.py` + - Needs: T011, T014 + - Wave: 3 + - Done when: `pytest tests/unit/modelling/test_fault_frame_analytic.py -k "inferred_centre or slip or extent_direction"` + passes. + - Context: `__init__`: `self._fault_buffer: float = 0.2`, + `self._force_mesh_geometry: bool = False`. `create_data_from_geometry` + stores both from its arguments first. `_on_geometry_change`: default + `force_mesh_geometry=None`; `None` args use the stored values. Inferred + centre: rows `coord == 0` and `val == 0`, columns X,Y,Z, `np.nanmean`; + a non-finite component -> 0.0 plus a warning naming the axis. Slip: + `s_perp = s - (s @ n) * n`; if `|s_perp| < 1e-12` raise + `ValueError("fault slip vector is parallel to the fault normal; supply a slip direction in the fault plane")`; + use the unit `s_perp` for `fault_slip_vector`, strike and depth rows. + `strike_unit = cross(n, s_perp)/|...|`. When `major_axis` is not None add + one row: X,Y,Z = centre, nx,ny,nz = `strike_unit / major_axis`, val NaN, + coord 2, w. Early return (no trace): call + `self.set_mesh_geometry(fault_buffer, None)` before `return`. + +- [x] T013 [US1] [S] Stop the centre projection and the buffer override in the fault factory. + - Files: `LoopStructural/modelling/core/_model_feature_factory.py`, + `LoopStructural/modelling/features/fault/_fault_segment.py` + - Needs: T006, T010, T012 + - Wave: 5 + - Done when: `pytest tests/unit/modelling/test_fault_frame_analytic.py tests/unit/modelling/test_fault_buffer_bounds.py` + passes. + - Context: In `create_and_add_fault`, delete the two lines + `if fault_center is not None and ~np.isnan(fault_center).any():` / + `fault_center = model.scale(fault_center, inplace=False)`. Delete the two + lines `if "force_mesh_geometry" not in kwargs:` / + `fault_frame_builder.set_mesh_geometry(kwargs.get("fault_buffer", 0.2), 0)`. + `create_data_from_geometry` already applies `fault_buffer` and + `force_mesh_geometry`. No signature change. An explicit `fault_center` + is a world point (FR-010). + Answer (plan.md "Answers", 2026-11-04): in `FaultSegment.fault_ellipsoid` + use `self.fault_centre[None, :]` with no `model.rescale`; the centre is a + world point. + +- [x] T014 [P] [US1] [M] Test: R5 buffer and forced-geometry bounds. + - Files: `tests/unit/modelling/test_fault_buffer_bounds.py` + - Needs: T002 + - Wave: 2 + - Done when: `pytest tests/unit/modelling/test_fault_buffer_bounds.py` + fails before T012 and T013, then passes after them. + - Context: Use `r5_fault_data`, `model_with_frame`, `rot_z`, `tp`. Box + O-(10,10,10) to O+(10,10,10); normal (1,0,0), slip (0,0,-1); axes major 12, + minor 4, intermediate 12; `force_mesh_geometry=True`; buffers b = 0, 0.2, + 0.5 through `create_and_add_fault`; FDI nelements 20000, PLI 1000; + explicit and inferred centre; rotations 0, 30, 90 deg. Requested world + box `O +/- (10+20*b)`; expected local box = enclosure of its 8 projected + corners. Tests: `test_buffer_bounds_match_request` (each support face is + outward of the expected bound by 0 to 1 cell (actual local spacing) plus + Tp; all generated rows inside the support); + `test_buffers_give_distinct_bounds`; + `test_property_change_keeps_buffer_and_forced_geometry` (set + `builder.fault_normal_vector = (1,0,0)` after creation; bounds unchanged). + +- [x] T015 [US1] [M] Align `test__fault_builder.py` with the corrected behaviour. + - Files: `tests/unit/modelling/test__fault_builder.py`, + `specs/005-coordinate-frame-correctness/verification.md` + - Needs: T001, T013 + - Wave: 6 + - Done when: `pytest tests/unit/modelling/test__fault_builder.py` passes. + - Context: Change an assertion only if it encodes a defect that S1 fixed + (for example an inferred centre with z = 0, a projected explicit centre, + or a default buffer that replaced the requested one). For each change add + one entry under `## Golden review (C-004)` in `verification.md`: test id, + old value, new value, defect, expected physical effect. Do not relax a + tolerance (SC-015). If no change is needed, write "No golden change in + S1." under that heading. + Answer (plan.md "Answers", 2026-11-04): also add one entry for the + T006 change to `packages/loop_common/tests/data/structured_support_reference.npz` + (27 values of `after_origin`/`after_maximum`; old values did not enclose + the requested bounds; FR-014). + +## Phase 5: Finish + +- [x] T016 [P] [Finish] [S] Add the S1 migration notice to `COMPAT.md`. + - Files: `COMPAT.md` + - Needs: T013 + - Wave: 6 + - Done when: the "Migration notices" table in `COMPAT.md` has the new row + with date 2026-11-04. + - Context: Add one row at the end of the table under + `## Migration notices (not breaking, no shim needed)`. Columns: Old | + New | Since. Old: "`create_and_add_fault(fault_center=...)` projected the + centre into local coordinates; `fault_buffer`/`force_mesh_geometry` were + replaced by a default 0.2 buffer; an inferred centre had z = 0; any 3x3 + matrix was accepted as a rotation". New: "`fault_center` is a world point; + the requested buffer and forced geometry are kept; an inferred centre keeps + its elevation; `set_local_transform` and `rotation_xy` raise + `RigidRotationError` for a non-rigid matrix. Extra observations added to + work around these defects are no longer necessary (spec 005)". Since: + 2026-11-04. + +- [x] T017 [Finish] [M] Record pre-fix evidence, gates, run time and rollback. + - Files: `specs/005-coordinate-frame-correctness/verification.md` + - Needs: T003, T004, T005, T006, T007, T008, T009, T010, T011, T012, T013, T014, T015, T016 + - Wave: 7 + - Done when: `ruff check LoopStructural packages && pytest` and + `pytest packages/loop_common/tests packages/loop_interpolation/tests` + pass, and each heading of `verification.md` has an entry. + - Context: Pre-fix evidence (FR-031): make a throw-away worktree at commit + `ebc7c699` (`git worktree add ebc7c699`), copy the 4 new test files + and `tests/fixtures/coordinate_frames.py` into it, run them with + `PYTHONPATH=;/packages/loop_common/src;/packages/loop_interpolation/src` + (the repo `.venv` is editable on the main checkout), record each failing + test id and value under `## Pre-fix failures`, then + `git worktree remove `. A finding that does not fail is "not + reproduced". Record pass results, the backend cells covered (Python FDI, + PLI/P1), the S1 test time (OP2), and the rollback check: on the + worktree at `ebc7c699`, `tests/unit/modelling/test__fault_builder.py` and + `packages/loop_common/tests/test_bounding_box.py` pass. + +## Wave table + +| Wave | Tasks | Tier [S] | Tier [M] | +|---|---|---|---| +| 1 | T001, T002, T003 | 1 | 2 | +| 2 | T004, T007, T011, T014 | 1 | 3 | +| 3 | T005, T006, T008, T012 | 1 | 3 | +| 4 | T009, T010 | 1 | 1 | +| 5 | T013 | 1 | 0 | +| 6 | T015, T016 | 1 | 1 | +| 7 | T017 | 0 | 1 | + +## Critical path + +T003 -> T004 -> T005 -> T010 -> T013 -> T015 -> T017 + +## Dependencies + +- Setup (T001, T002) and the foundation test T003 start at once. +- Foundation: T003 -> T004 -> T005, T006. +- US2: T007 -> T008 -> T009, T010 (T010 also needs T005). +- US1: T011, T014 -> T012 -> T013 (T013 also needs T006 and T010, because + the fault frame uses the repaired reference sampling and resize). +- Finish: T015 and T016 after T013; T017 last. +- Stages S2 to S6 have 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 new file mode 100644 index 000000000..5bfc0e914 --- /dev/null +++ b/specs/005-coordinate-frame-correctness/verification.md @@ -0,0 +1,137 @@ +# Verification: 005-coordinate-frame-correctness + +This file records the evidence for FR-029, FR-031 and FR-032 of spec 005. + +## Pre-fix failures + +Method. A throw-away worktree at `ebc7c699` (the commit before the spec) held +only the 4 new test files and `tests/fixtures/coordinate_frames.py`. The +`tests/__init__.py` and `tests/fixtures/__init__.py` files already exist at +`ebc7c699`, so no `__init__.py` was copied. The run used +`PYTHONPATH=;/packages/loop_common/src;/packages/loop_interpolation/src`. +Before each run, `python -c "import loop_common, LoopStructural"` showed that +both imports came from the throw-away worktree. The worktree is removed +(`git worktree remove --force`, then `git worktree prune`). + +| Test file | Result at `ebc7c699` | Time | +|---|---|---| +| `packages/loop_common/tests/test_rigid_transform_contract.py` | Collection error: `ImportError: cannot import name 'RigidRotationError' from 'loop_common.geometry'`. No test runs. | 2 s | +| `tests/unit/modelling/test_reference_constraint_frames.py` | 278 failed, 34 passed | 7 s | +| `tests/unit/modelling/test_fault_frame_analytic.py` | 53 failed, 0 passed | 88 s | +| `tests/unit/modelling/test_fault_buffer_bounds.py` | 84 failed, 0 passed | 12 s | + +### Rigid-transform contract (no R1 fixture; FR-003, FR-013, FR-014, FR-021; SC-001, SC-011, SC-012) +The file cannot import pre-fix, because `RigidRotationError` does not exist. +All 42 contract checks are "not runnable" pre-fix. None is "reproduced" or +"not reproduced" by a test result. The missing class is itself the evidence. + +### Reference constraints (R3; FR-003 to FR-007, FR-009, FR-013, FR-024; SC-007, SC-008) +| Test function | Failed / passed | Representative failing id and value | +|---|---|---| +| `test_world_local_conversion_methods` | 24 / 0 | `[frame_case0-FDI]`: `AttributeError: 'FiniteDifferenceInterpolator' object has no attribute 'world_to_local_points'` (same for `P1Interpolator`) | +| `test_reference_gradient_samples_use_world_points_and_local_vectors` | 66 / 6 | first mismatch at `[0, 0]`: -3.807 (actual) vs -2.770 (expected), atol 1e-10 | +| `test_orthogonal_rows_affine_residual` | 48 / 24 | `[frame_case1-None-FDI]`: 210 of 210 residuals about -1.687 vs 0.0 | +| `test_inequality_feature_bounds_world` | 44 / 4 | `[frame_case1-True-FDI]`: bound -4.9999 vs -4.0176 (336 of 336 differ) | +| `test_equality_region_selects_world_nodes` | 24 / 0 | `'reference_equality'` is not in the constraint dictionary (rows are under `'undefined'`) | +| `test_restricted_reference_exclusion_counts` | 24 / 0 | `'restricted_equality'` is not in the constraint dictionary | +| `test_set_interpolation_geometry_encloses_rotated_box` | 24 / 0 | `[frame_case0-FDI]`: grid maximum [6.4, 8.33, 10.29] does not enclose the box corner [0, 10, 12]; also `ValueError: Cannot resize the grid` in some cells | +| `test_support_expansion_contains_restored_points` | 24 / 0 | `[frame_case0-FDI]`: grid maximum [27.2, 8.67, 13.71] does not enclose the restored points | + +### Fault frame analytic (R2, R4; FR-001, FR-002, FR-010, FR-011, FR-026; SC-003, SC-013) +All 5 single tests fail. The 48 planar cells all fail (12 for each of +FDI/PLI x centre mode). + +| Test function | Failed / passed | Representative failing id and value | +|---|---|---| +| `test_inferred_centre_keeps_elevation` | 1 / 0 | inferred centre [100, 200, 0] vs expected [100, 200, 300] | +| `test_explicit_centre_rows_not_projected` | 1 / 0 | no explicit centre row equals the centre [100, 200, 300]; the rows were projected | +| `test_slip_vector_made_orthogonal` | 1 / 0 | slip [0.7071, 0, -0.7071] vs expected [0, 0, -1] | +| `test_parallel_slip_raises` | 1 / 0 | `DID NOT RAISE ValueError` | +| `test_extent_direction_row_at_centre` | 1 / 0 | the count of rows at the centre is 0, expected 1 | +| `test_planar_fault_frame_matches_analytic` | 48 / 0 | `[FDI-True-0.0-anchor0]`: maximum field error 0.830 (limit 1e-6). Error range over cells: 1.7e-3 to 0.84. 12 cells show exactly 0.45. | + +### Fault buffer bounds (R5; FR-012, FR-013; SC-005) +| Test function | Failed / passed | Representative failing id and value | +|---|---|---| +| `test_buffer_bounds_match_request` | 36 / 0 | `[0.0-True-0.0-FDI-20000]`: lower bound 0 vs the allowed -4 (bounds do not cover the buffer) | +| `test_buffers_give_distinct_bounds` | 12 / 0 | `[True-0.0-FDI-20000]`: same bound failure | +| `test_property_change_keeps_buffer_and_forced_geometry` | 36 / 0 | `ValueError: Cannot resize the grid` ("proposed number of steps is [-2 13 13]") and bound assertion failures | + +### Findings that are "not reproduced" +None in the three runnable files: each test function has at least one +failing cell. The partial passes (for example 24 of 72 in +`test_orthogonal_rows_affine_residual`) are cells where the case does not +exercise the defect. + +## Post-fix results + +Spec worktree, HEAD `10fd4dea`, worktree `PYTHONPATH`. A check of +`loop_common.__file__` and `LoopStructural.__file__` showed that both imports +come from `C:\Users\lachl\dev\LoopStructural-005`. + +| Command | Result | +|---|---| +| `python -m pytest packages/loop_common/tests packages/loop_interpolation/tests -p no:cacheprovider -q` | 912 passed, 47 skipped, 2 xfailed, 7 warnings in 142 s | +| The 4 S1 test files in one run | 491 passed, 7 warnings in 291 s | +| `test_rigid_transform_contract.py` | 42 passed | +| `test_reference_constraint_frames.py` | 312 passed | +| `test_fault_frame_analytic.py` | 53 passed | +| `test_fault_buffer_bounds.py` | 84 passed | +| `ruff check LoopStructural packages` (orchestrator) | All checks passed | +| Root `pytest` (orchestrator, full) | 1204 passed, 5 skipped, 21 warnings in 1082.66 s. Baseline at `3881ef32`: 755 passed, 5 skipped, 21 warnings in 1205.98 s. | + +Test-setup corrections made during the run (facts for the record). No +tolerance was relaxed. +- T003: the round-trip test passes `dimensions=dimensions`. +- T011: the slip-row expectation is unit slip / intermediate_axis (spec R2 table). +- T007: the affine residual computes the exact field without large-coordinate + cancellation. The atol of 1e-10 is kept. +- T011: the helper passes the SC-003 CG settings through `solver_kwargs`. + Top-level values were ignored, so the solver used atol 2e-9. + +## Golden review (C-004) + +No golden change in S1 for `tests/unit/modelling/test__fault_builder.py`. +No assertion in that file encodes an S1 defect. The file checks the +`update_geometry` bounds, and the fault normal and slip vectors. It does not +check a fault centre or a default buffer. The commented centre assertions +stay unchanged. + +One golden change is made in S1 (task T006, commit d2909aa5). + +| Item | Old value | New value | Defect | Expected physical effect | +|---|---|---|---|---| +| `packages/loop_common/tests/data/structured_support_reference.npz`: 27 values `//after_origin.nsteps`, `after_maximum.nsteps`, `after_maximum.maximum` for cls StructuredGrid, TetMesh, P2TetMesh and grid unit, aniso, small | unit: origin nsteps [4 4 4], maximum nsteps [6 6 6], maximum [6.142857 x3]. aniso: [3 6 4], [5 7 5], [14.6 6.909091 7.857143]. small: [5 5 5], [8 8 8], [5.0 x3] | unit: [5 5 5], [8 8 8], [9.0 x3]. aniso: [4 7 5], [7 9 7], [17.4 10.545455 11.285714]. small: [6 6 6], [10 10 10], [6.714286 x3] | The origin and maximum setters stored a cell count as the node count. The support lost one cell and did not reach the requested maximum (FR-014). Old values did not enclose the requested bounds (for example unit: requested maximum [7.714286 x3]). | The support now encloses the requested bounds, within one cell, outward. | + +The other keys of the file keep their values from before S1 (the 33 keys that `--update` re-recorded with round-off changes were restored after the spec review). No tolerance was relaxed. + + +## Backend matrix + +S1 tests use Python backends only. + +| Requirement | Backends | Cells | +|---|---|---| +| R2 fault frame analytic | Python FDI, Python PLI (P1) | 4 anchors x 3 rotations x 2 centre modes x 2 backends = 48 planar cells, plus 5 single tests | +| R3 reference constraints | FDI, PLI (P1) | reference file, sample rules None / centroid / gauss4, all frame cases (312 cells) | +| R5 fault buffers | FDI, PLI | buffers 0, 0.2, 0.5 x 3 rotations x 2 centre modes (84 cells) | + +No native, browser or RBF backend is in S1 scope. FR-030 states that such +backends are not a prerequisite. + +## Rollback checks + +On a throw-away worktree at `ebc7c699` (imports checked as above): +`python -m pytest tests/unit/modelling/test__fault_builder.py packages/loop_common/tests/test_bounding_box.py -p no:cacheprovider -q` +gives 15 passed in 0.5 s. The two files pass on the pre-spec code, so the +rollback target is sound. + +## Test run time (OP2) + +| Item | Time | +|---|---| +| All 4 S1 files, one run | 291 s test time (4 min 57 s wall) | +| Slowest file: `test_fault_frame_analytic.py` | 264 s test time (4 min 31 s wall) | +| `test_fault_buffer_bounds.py` | 20 s (26 s wall) | +| `test_reference_constraint_frames.py` | 5.6 s (12 s wall) | +| `test_rigid_transform_contract.py` | 0.1 s (4 s wall) | diff --git a/specs/README.md b/specs/README.md index 0288aeb37..d2c592603 100644 --- a/specs/README.md +++ b/specs/README.md @@ -1,6 +1,6 @@ # Specification index -Updated: 2026-11-03. This index records document and workflow state. It does +Updated: 2026-11-04. This index records document and workflow state. It does not certify that an implementation or a numerical acceptance gate passed. | Spec | State | Next step or existing open point | @@ -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) | Cross-spec documentation prerequisite reconciled | Continue the 005 planning workflow after confirming SC-019. Its source header still records the earlier cascade block; this cascade leaves the source spec unchanged. | +| [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)); S2 to S6 not planned | Review and merge the S1 branch. OP1: plan S2 to S6 in later planner runs. 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/fixtures/coordinate_frames.py b/tests/fixtures/coordinate_frames.py new file mode 100644 index 000000000..3bf32fa74 --- /dev/null +++ b/tests/fixtures/coordinate_frames.py @@ -0,0 +1,176 @@ +"""Shared analytic data and reference functions for coordinate frame tests.""" + +from __future__ import annotations + +from collections.abc import Callable + +import numpy as np +import pandas as pd + +from LoopStructural import GeologicalModel + +ANCHORS = np.array( + [ + (0.0, 0.0, 0.0), + (100.0, 200.0, 300.0), + (-100.0, -200.0, -300.0), + (500000.0, 6500000.0, 100.0), + ] +) +ROTATIONS_DEG = (0.0, 30.0, 90.0) + + +def rot_z(deg: float) -> np.ndarray: + """Return a rotation matrix for a rotation around the z axis.""" + angle = np.deg2rad(deg) + cosine, sine = np.cos(angle), np.sin(angle) + return np.array([[cosine, -sine, 0.0], [sine, cosine, 0.0], [0.0, 0.0, 1.0]]) + + +def tp(length: float, width: float) -> float: + """Return a scale-aware floating point tolerance.""" + return 1e-10 * length + 32 * np.finfo(float).eps * max(width, length) + + +def model_with_frame( + origin: np.ndarray, maximum: np.ndarray, rotation_deg: float +) -> GeologicalModel: + """Create a model with a rotated local bounding-box frame.""" + model = GeologicalModel(origin, maximum) + model.bounding_box.set_local_transform(local_origin=origin, rotation_matrix=rot_z(rotation_deg)) + return model + + +def _fault_data(anchor: np.ndarray, offsets: tuple[int, ...]) -> pd.DataFrame: + y_values, z_values = np.meshgrid(offsets, offsets, indexing="ij") + local = np.column_stack((np.zeros(y_values.size), y_values.ravel(), z_values.ravel())) + points = anchor + local + count = points.shape[0] + return pd.DataFrame( + { + "X": points[:, 0], + "Y": points[:, 1], + "Z": points[:, 2], + "val": np.zeros(count), + "coord": np.zeros(count), + "feature_name": ["fault"] * count, + "w": np.ones(count), + } + ) + + +def r2_fault_data(anchor: np.ndarray) -> pd.DataFrame: + """Return fault observations for requirement R2.""" + return _fault_data(anchor, (-8, -4, 0, 4, 8)) + + +def r5_fault_data(anchor: np.ndarray) -> pd.DataFrame: + """Return fault observations for requirement R5.""" + return _fault_data(anchor, (-4, -2, 0, 2, 4)) + + +def r2_probes(anchor: np.ndarray) -> np.ndarray: + """Return the 19-cubed probe points around an anchor.""" + offsets = np.arange(-9, 10, dtype=float) + grid = np.stack(np.meshgrid(offsets, offsets, offsets, indexing="ij"), axis=-1) + return anchor + grid.reshape(-1, 3) + + +def r2_frame(points: np.ndarray, anchor: np.ndarray) -> tuple[np.ndarray, np.ndarray]: + """Evaluate the R2 frame values and gradients at world points.""" + local = np.asarray(points) - anchor + x, y, z = local.T + values = np.column_stack((x / 20, -z / 20, y / 20)) + gradients = np.broadcast_to( + np.array([[1 / 20, 0, 0], [0, 0, -1 / 20], [0, 1 / 20, 0]]), + (len(local), 3, 3), + ).copy() + return values, gradients + + +def r3_value(points: np.ndarray, anchor: np.ndarray) -> np.ndarray: + """Evaluate the R3 analytic scalar field.""" + local = np.asarray(points) - anchor + x, y, z = local.T + return x + 0.05 * y**2 + z + + +def r3_gradient(points: np.ndarray, anchor: np.ndarray) -> np.ndarray: + """Evaluate the R3 analytic scalar field gradient.""" + local = np.asarray(points) - anchor + gradient = np.zeros((len(local), 3), dtype=float) + gradient[:, 0] = 1.0 + gradient[:, 1] = 0.1 * local[:, 1] + gradient[:, 2] = 1.0 + return gradient + + +def r4_frame( + points: np.ndarray, anchor: np.ndarray +) -> tuple[np.ndarray, tuple[np.ndarray, np.ndarray, np.ndarray]]: + """Evaluate the R4 cylindrical frame values and gradients.""" + local = np.asarray(points) - anchor + x, y, z = local.T + radius_origin = 40.0 + world_x = radius_origin + x + radius = np.sqrt(world_x**2 + y**2) + values = np.column_stack( + ( + (radius - radius_origin) / 20, + -z / 20, + radius_origin * np.arctan2(y, world_x) / 20, + ) + ) + gradients = ( + np.column_stack((world_x / (20 * radius), y / (20 * radius), np.zeros(len(local)))), + np.broadcast_to([0.0, 0.0, -1 / 20], (len(local), 3)).copy(), + np.column_stack( + ( + -radius_origin * y / (20 * radius**2), + radius_origin * world_x / (20 * radius**2), + np.zeros(len(local)), + ) + ), + ) + return values, gradients + + +class AnalyticReference: + """Evaluate an analytic scalar field and gradient within an optional box.""" + + def __init__( + self, + value_fn: Callable[[np.ndarray], np.ndarray], + gradient_fn: Callable[[np.ndarray], np.ndarray], + name: str, + support_box: tuple[np.ndarray, np.ndarray] | None = None, + ) -> None: + self.value_fn = value_fn + self.gradient_fn = gradient_fn + self.name = name + self.support_box = support_box + + def _supported(self, points: np.ndarray) -> np.ndarray: + if self.support_box is None: + return np.ones(len(points), dtype=bool) + lower, upper = self.support_box + return np.all((points >= lower) & (points <= upper), axis=1) + + def evaluate_value(self, points: np.ndarray) -> np.ndarray: + """Evaluate values, returning NaN outside the closed support box.""" + points = np.asarray(points) + values = np.asarray(self.value_fn(points), dtype=float).copy() + values[~self._supported(points)] = np.nan + return values + + def evaluate_gradient(self, points: np.ndarray, ignore_regions: bool = True) -> np.ndarray: + """Evaluate gradients, returning NaN outside the closed support box.""" + del ignore_regions + points = np.asarray(points) + gradients = np.asarray(self.gradient_fn(points), dtype=float).copy() + gradients[~self._supported(points)] = np.nan + return gradients + + def __call__(self, points: np.ndarray) -> np.ndarray: + """Evaluate the reference values.""" + return self.evaluate_value(points) diff --git a/tests/unit/modelling/test_fault_buffer_bounds.py b/tests/unit/modelling/test_fault_buffer_bounds.py new file mode 100644 index 000000000..28733e8a0 --- /dev/null +++ b/tests/unit/modelling/test_fault_buffer_bounds.py @@ -0,0 +1,115 @@ +"""Check buffered fault supports in rotated model coordinate frames.""" + +from itertools import product + +import numpy as np +import pytest + +from tests.fixtures.coordinate_frames import ROTATIONS_DEG, model_with_frame, r5_fault_data, tp + +ORIGIN = np.array([500000.0, 6500000.0, 100.0]) +HALF_SIZE = np.full(3, 10.0) +BUFFERS = (0.0, 0.2, 0.5) +ENGINES = (("FDI", 20_000), ("PLI", 1_000)) + + +def _make_fault(buffer, rotation_deg, explicit_center, engine, nelements): + model = model_with_frame(ORIGIN - HALF_SIZE, ORIGIN + HALF_SIZE, rotation_deg) + model.data = r5_fault_data(ORIGIN) + fault = model.create_and_add_fault( + "fault", + 1.0, + data=model.data, + interpolatortype=engine, + nelements=nelements, + fault_normal_vector=np.array([1.0, 0.0, 0.0]), + fault_slip_vector=np.array([0.0, 0.0, -1.0]), + fault_center=ORIGIN if explicit_center else None, + major_axis=12, + minor_axis=4, + intermediate_axis=12, + force_mesh_geometry=True, + fault_buffer=buffer, + ) + return model, fault.builder + + +def _expected_bounds(model, buffer): + radius = 10.0 + 20.0 * buffer + world_corners = np.array( + list(product(*[(ORIGIN[i] - radius, ORIGIN[i] + radius) for i in range(3)])) + ) + local_corners = model.bounding_box.project(world_corners) + return np.min(local_corners, axis=0), np.max(local_corners, axis=0), world_corners + + +def _supports(builder): + return [coordinate_builder.interpolator.support for coordinate_builder in builder.builders] + + +def _assert_support_matches(model, builder, buffer): + expected_origin, expected_maximum, world_corners = _expected_bounds(model, buffer) + length = np.max(expected_maximum - expected_origin) + world_magnitude = np.max(np.abs(world_corners)) + for coordinate_builder, support in zip(builder.builders, _supports(builder)): + spacing = np.asarray(support.step_vector) + tolerance = tp(length, world_magnitude) + actual_origin = np.asarray(support.origin) + actual_maximum = np.asarray(support.maximum) + assert np.all(expected_origin - actual_origin >= -tolerance) + assert np.all(expected_origin - actual_origin <= spacing + tolerance) + assert np.all(actual_maximum - expected_maximum >= -tolerance) + assert np.all(actual_maximum - expected_maximum <= spacing + tolerance) + + rows = coordinate_builder.data[["X", "Y", "Z"]].to_numpy(dtype=float) + rows = rows[np.isfinite(rows).all(axis=1)] + assert rows.size > 0 + local_rows = coordinate_builder.interpolator.bounding_box.project(rows) + assert np.all(local_rows >= actual_origin - tolerance) + assert np.all(local_rows <= actual_maximum + tolerance) + + +@pytest.mark.parametrize("engine,nelements", ENGINES) +@pytest.mark.parametrize("rotation_deg", ROTATIONS_DEG) +@pytest.mark.parametrize("explicit_center", [True, False]) +@pytest.mark.parametrize("buffer", BUFFERS) +def test_buffer_bounds_match_request(buffer, explicit_center, rotation_deg, engine, nelements): + model, builder = _make_fault(buffer, rotation_deg, explicit_center, engine, nelements) + _assert_support_matches(model, builder, buffer) + + +@pytest.mark.parametrize("engine,nelements", ENGINES) +@pytest.mark.parametrize("rotation_deg", ROTATIONS_DEG) +@pytest.mark.parametrize("explicit_center", [True, False]) +def test_buffers_give_distinct_bounds(explicit_center, rotation_deg, engine, nelements): + bounds = [] + for buffer in BUFFERS: + model, builder = _make_fault(buffer, rotation_deg, explicit_center, engine, nelements) + supports = _supports(builder) + bounds.append([(support.origin.copy(), support.maximum.copy()) for support in supports]) + _assert_support_matches(model, builder, buffer) + + for index in range(3): + assert not np.allclose(bounds[0][index][0], bounds[1][index][0]) + assert not np.allclose(bounds[1][index][0], bounds[2][index][0]) + assert not np.allclose(bounds[0][index][1], bounds[1][index][1]) + assert not np.allclose(bounds[1][index][1], bounds[2][index][1]) + + +@pytest.mark.parametrize("engine,nelements", ENGINES) +@pytest.mark.parametrize("rotation_deg", ROTATIONS_DEG) +@pytest.mark.parametrize("explicit_center", [True, False]) +@pytest.mark.parametrize("buffer", BUFFERS) +def test_property_change_keeps_buffer_and_forced_geometry( + buffer, explicit_center, rotation_deg, engine, nelements +): + model, builder = _make_fault(buffer, rotation_deg, explicit_center, engine, nelements) + before = [(support.origin.copy(), support.maximum.copy()) for support in _supports(builder)] + + builder.fault_normal_vector = np.array([1.0, 0.0, 0.0]) + + after = _supports(builder) + for (origin_before, maximum_before), support in zip(before, after): + np.testing.assert_array_equal(support.origin, origin_before) + np.testing.assert_array_equal(support.maximum, maximum_before) + _assert_support_matches(model, builder, buffer) diff --git a/tests/unit/modelling/test_fault_frame_analytic.py b/tests/unit/modelling/test_fault_frame_analytic.py new file mode 100644 index 000000000..6a1d432bf --- /dev/null +++ b/tests/unit/modelling/test_fault_frame_analytic.py @@ -0,0 +1,229 @@ +"""Check analytically defined fault frames in model coordinate frames.""" + +import numpy as np +import pytest + +from tests.fixtures.coordinate_frames import ( + ANCHORS, + ROTATIONS_DEG, + model_with_frame, + r2_fault_data, + r2_frame, + r2_probes, + r4_frame, +) + +HALF_SIZE = np.full(3, 10.0) +NORMAL = np.array([1.0, 0.0, 0.0]) +SLIP = np.array([0.0, 0.0, -1.0]) + + +def _fault(anchor, rotation_deg, interpolatortype, explicit_center=True, slip=SLIP): + model = model_with_frame(anchor - HALF_SIZE, anchor + HALF_SIZE, rotation_deg) + data = r2_fault_data(anchor) + model.data = data + fault = model.create_and_add_fault( + "fault", + 0, + data=data, + interpolatortype=interpolatortype, + fault_center=anchor if explicit_center else None, + fault_normal_vector=NORMAL, + fault_slip_vector=slip, + major_axis=20, + minor_axis=10, + intermediate_axis=20, + fault_buffer=0.2, + nelements=20_000 if interpolatortype == "FDI" else 1_000, + fault_trace_anisotropy=0, + fault_dip_anisotropy=0, + ) + # These are the CG solve settings of SC-003. The builder ignores + # top-level rtol, atol and maxiter, so they go in solver_kwargs. + fault.builder.build( + solver="cg", + solver_kwargs={"rtol": 1e-10, "atol": 1e-12, "maxiter": 20_000}, + w1=1, + w2=1, + w3=1, + step=10, + ) + return fault + + +def test_explicit_centre_rows_not_projected(): + anchor = ANCHORS[1] + fault = _fault(anchor, 30, "FDI") + for builder in fault.builder.builders: + generated = builder.data + rows = generated.loc[generated["val"] == 0, ["X", "Y", "Z"]].to_numpy() + assert np.any(np.all(rows == anchor, axis=1)) + + +def test_inferred_centre_keeps_elevation(): + anchor = ANCHORS[1] + fault = _fault(anchor, 30, "FDI", explicit_center=False) + np.testing.assert_allclose(fault.builder.fault_centre, anchor, rtol=0, atol=1e-12) + assert fault.builder.fault_centre[2] == anchor[2] + + +def test_slip_vector_made_orthogonal(): + fault = _fault(ANCHORS[0], 0, "FDI", slip=np.array([1.0, 0.0, -1.0])) + normal = fault.builder.fault_normal_vector + slip = fault.builder.fault_slip_vector + np.testing.assert_allclose(slip, [0, 0, -1], rtol=0, atol=1e-12) + assert abs(np.dot(normal, slip)) <= 1e-12 + coordinate1 = fault.builder.builders[1].data + gradients = coordinate1[["nx", "ny", "nz"]].to_numpy() + # The coordinate 1 gradient is the unit in-plane slip divided by the intermediate axis. + expected = np.array([0.0, 0.0, -1.0]) / fault.builder.fault_intermediate_axis + assert np.any(np.all(np.isclose(gradients, expected, rtol=0, atol=1e-12), axis=1)) + + +def test_parallel_slip_raises(): + model = model_with_frame(ANCHORS[0] - HALF_SIZE, ANCHORS[0] + HALF_SIZE, 0) + data = r2_fault_data(ANCHORS[0]) + with pytest.raises(ValueError): + model.create_and_add_fault( + "fault", + 0, + data=data, + fault_center=ANCHORS[0], + fault_normal_vector=NORMAL, + fault_slip_vector=NORMAL, + major_axis=20, + minor_axis=10, + intermediate_axis=20, + ) + + +def test_extent_direction_row_at_centre(): + anchor = ANCHORS[0] + fault = _fault(anchor, 0, "FDI", explicit_center=False) + coordinate2 = fault.builder.builders[2].data + directions = coordinate2[["nx", "ny", "nz"]].to_numpy() + finite = np.all(np.isfinite(directions), axis=1) + assert np.count_nonzero(finite) == 1 + assert np.all(coordinate2.loc[finite, ["X", "Y", "Z"]].to_numpy()[0] == anchor) + row = coordinate2.loc[finite] + np.testing.assert_allclose( + row[["nx", "ny", "nz"]].to_numpy()[0], [0, 1 / 20, 0], rtol=0, atol=1e-10 + ) + _, gradients = r4_frame(anchor[None, :], anchor) + np.testing.assert_allclose( + row[["nx", "ny", "nz"]].to_numpy()[0], gradients[2][0], rtol=0, atol=1e-10 + ) + points = anchor + np.array([[0, -6, 0], [0, 6, 0]]) + directions = r4_frame(points, anchor)[1][2] + cosine = ( + np.dot(directions[0], directions[1]) + / np.linalg.norm(directions[0]) + / np.linalg.norm(directions[1]) + ) + assert np.rad2deg(np.arccos(np.clip(cosine, -1, 1))) >= 15 + + +@pytest.mark.parametrize("anchor", ANCHORS.tolist()) +@pytest.mark.parametrize("rotation_deg", ROTATIONS_DEG) +@pytest.mark.parametrize("explicit_center", [True, False]) +@pytest.mark.parametrize("interpolatortype", ["FDI", "PLI"]) +def test_planar_fault_frame_matches_analytic( + anchor, rotation_deg, explicit_center, interpolatortype +): + anchor = np.asarray(anchor) + fault = _fault(anchor, rotation_deg, interpolatortype, explicit_center) + points = r2_probes(anchor) + expected, expected_gradients = r2_frame(points, anchor) + actual = np.column_stack([fault[i].evaluate_value(points) for i in range(3)]) + actual_gradients = np.stack([fault[i].evaluate_gradient(points) for i in range(3)], axis=1) + assert np.all(np.isfinite(actual)) + assert np.all(np.isfinite(actual_gradients)) + assert np.max(np.abs(actual - expected)) <= 1e-6 + expected_magnitude = np.linalg.norm(expected_gradients, axis=2) + actual_magnitude = np.linalg.norm(actual_gradients, axis=2) + assert ( + np.max(np.abs(actual_magnitude - expected_magnitude)) / np.max(expected_magnitude) <= 1e-4 + ) + dot = np.sum(actual_gradients * expected_gradients, axis=2) + denom = np.linalg.norm(actual_gradients, axis=2) * expected_magnitude + angles = np.rad2deg(np.arccos(np.clip(dot / denom, -1, 1))) + assert np.all(np.isfinite(angles)) + assert np.max(angles) <= 0.01 + + +def test_caller_normal_not_mutated(): + # The builder must copy the normal vector of the caller. + anchor = ANCHORS[1] + model = model_with_frame(anchor - HALF_SIZE, anchor + HALF_SIZE, 30) + data = r2_fault_data(anchor) + model.data = data + normal = np.array([2.0, 0.0, 0.0]) + slip = np.array([0.0, 0.0, -2.0]) + model.create_and_add_fault( + "fault", + 0, + data=data, + interpolatortype="FDI", + fault_center=anchor, + fault_normal_vector=normal, + fault_slip_vector=slip, + major_axis=20, + minor_axis=10, + intermediate_axis=20, + fault_buffer=0.2, + nelements=20_000, + fault_trace_anisotropy=0, + fault_dip_anisotropy=0, + ) + np.testing.assert_array_equal(normal, [2.0, 0.0, 0.0]) + np.testing.assert_array_equal(slip, [0.0, 0.0, -2.0]) + + +def test_fault_ellipsoid_centre_is_world(monkeypatch): + # fault_centre is a world point, so the ellipsoid must not convert it again. + import sys + import types + + received = {} + + class FakePolyData: + def __init__(self, pts): + received["pts"] = np.array(pts, dtype=float) + self._data = {} + + def __setitem__(self, key, value): + self._data[key] = value + + def __getitem__(self, key): + return self._data[key] + + def glyph(self, **kwargs): + return self + + fake = types.ModuleType("pyvista") + fake.PolyData = FakePolyData + fake.ParametricEllipsoid = lambda *args: object() + monkeypatch.setitem(sys.modules, "pyvista", fake) + + anchor = ANCHORS[1] + model = model_with_frame(anchor - HALF_SIZE, anchor + HALF_SIZE, 30) + data = r2_fault_data(anchor) + model.data = data + fault = model.create_and_add_fault( + "fault", + 0, + data=data, + interpolatortype="FDI", + fault_center=anchor, + fault_normal_vector=NORMAL, + fault_slip_vector=SLIP, + major_axis=20, + minor_axis=10, + intermediate_axis=20, + fault_buffer=0.2, + nelements=20_000, + fault_trace_anisotropy=0, + fault_dip_anisotropy=0, + ) + fault.fault_ellipsoid() + np.testing.assert_allclose(received["pts"][0], anchor, atol=1e-9) diff --git a/tests/unit/modelling/test_reference_constraint_frames.py b/tests/unit/modelling/test_reference_constraint_frames.py new file mode 100644 index 000000000..5c33694f8 --- /dev/null +++ b/tests/unit/modelling/test_reference_constraint_frames.py @@ -0,0 +1,246 @@ +"""Reference frame checks for discrete interpolation constraints.""" + +import itertools + +import numpy as np +import pytest +from loop_common.supports import StructuredGrid, TetMesh + +from LoopStructural.modelling.features.builders import GeologicalFeatureBuilder +from tests.fixtures.coordinate_frames import ( + ANCHORS, + ROTATIONS_DEG, + AnalyticReference, + model_with_frame, + r3_gradient, + r3_value, + tp, +) + + +@pytest.fixture(params=list(itertools.product(range(len(ANCHORS)), ROTATIONS_DEG))) +def frame_case(request): + anchor, rotation = ANCHORS[request.param[0]], request.param[1] + model = model_with_frame(anchor, anchor + (8, 10, 12), rotation) + return anchor, rotation, model + + +def _builder(engine, model): + builder = GeologicalFeatureBuilder( + interpolatortype=engine, + bounding_box=model.bounding_box, + nelements=210, + name="r3", + model=model, + buffer=0, + ) + corners = np.array( + list(itertools.product(*zip(model.bounding_box.origin, model.bounding_box.maximum))) + ) + local = model.bounding_box.project(corners) + lo, hi = local.min(axis=0), local.max(axis=0) + support_type = StructuredGrid if engine == "FDI" else TetMesh + interp = type(builder.interpolator)( + support_type(origin=lo, nsteps_cells=[5, 6, 7], step_vector=(hi - lo) / [5, 6, 7]) + ) + interp.bounding_box = model.bounding_box + builder._interpolator = interp + return builder + + +def _field(anchor): + return AnalyticReference( + lambda p: r3_value(p, anchor), lambda p: r3_gradient(p, anchor), "reference" + ) + + +@pytest.mark.parametrize("engine", ["FDI", "PLI"]) +def test_world_local_conversion_methods(frame_case, engine): + anchor, _, model = frame_case + interp = _builder(engine, model).interpolator + points = anchor + np.array([[1.2, 2.3, 3.4], [4.1, 5.2, 6.3]]) + vectors = np.array([[1.0, 2.0, 3.0], [-2.0, 1.0, 4.0]]) + for to_local, to_world, values in ( + (interp.world_to_local_points, interp.local_to_world_points, points), + (interp.world_to_local_vectors, interp.local_to_world_vectors, vectors), + ): + np.testing.assert_allclose(to_world(to_local(values)), values, atol=tp(12, 12), rtol=0) + + +@pytest.mark.parametrize("engine", ["FDI", "PLI"]) +def test_equality_region_selects_world_nodes(frame_case, engine): + anchor, _, model = frame_case + builder = _builder(engine, model) + nodes = builder.interpolator.support.nodes + world = model.bounding_box.reproject(nodes) + expected = np.flatnonzero(world[:, 0] > anchor[0] + 4) + field = _field(anchor) + builder.add_equality_constraints(field, lambda p: p[:, 0] > anchor[0] + 4) + builder.install_equality_constraints() + entries = builder.interpolator.equal_constraints + actual = np.concatenate([entry["col"] for entry in entries.values()]) + np.testing.assert_array_equal(np.sort(actual), expected) + np.testing.assert_allclose( + np.concatenate([entry["B"] for entry in entries.values()]), + field.evaluate_value(world[expected]), + rtol=0, + atol=tp(50, 50), + ) + assert "reference_equality" in entries + + +@pytest.mark.parametrize("engine", ["FDI", "PLI"]) +@pytest.mark.parametrize("sample_rule", [None, "centroid", "gauss4"]) +def test_orthogonal_rows_affine_residual(frame_case, engine, sample_rule): + anchor, _, model = frame_case + builder = _builder(engine, model) + for name, world_direction, target in ( + ("x", np.array([1.0, 0.0, 0.0]), 0.0), + ("y", np.array([0.0, 1.0, 0.0]), 1.0), + ): + reference = AnalyticReference( + lambda p, offset=anchor[1]: p[:, 1] - offset, + lambda p, direction=world_direction: np.broadcast_to(direction, (len(p), 3)).copy(), + name, + ) + builder.add_orthogonal_feature(reference, B=target, sample_rule=sample_rule) + builder.install_gradient_constraint() + support = builder.interpolator.support + bounding_box = model.bounding_box + # Use the rigid structure to avoid cancellation at large world coordinates. + exact = bounding_box.reproject_vectors(support.nodes)[:, 1] + ( + bounding_box.local_origin[1] - anchor[1] + ) + for name in ("x_orthogonal", "y_orthogonal"): + constraint = builder.interpolator.constraints[name] + residual = constraint["matrix"] @ exact - constraint["b"] + assert len(residual) > 0 + np.testing.assert_allclose(residual, 0.0, rtol=0, atol=1e-10) + + +@pytest.mark.parametrize("engine", ["FDI", "PLI"]) +@pytest.mark.parametrize("sample_rule", [None, "centroid", "gauss4"]) +def test_reference_gradient_samples_use_world_points_and_local_vectors( + frame_case, engine, sample_rule +): + anchor, _, model = frame_case + builder = _builder(engine, model) + reference = _field(anchor) + seen = {} + original_eval = reference.evaluate_gradient + + def record_gradient(points, ignore_regions=True): + seen["points"] = np.asarray(points).copy() + return original_eval(points, ignore_regions=ignore_regions) + + reference.evaluate_gradient = record_gradient + original_add = builder.interpolator.add_gradient_orthogonal_constraints + + def record_add(points, vector, **kwargs): + seen["vectors"] = np.asarray(vector).copy() + return original_add(points, vector, **kwargs) + + builder.interpolator.add_gradient_orthogonal_constraints = record_add + builder.add_orthogonal_feature(reference, sample_rule=sample_rule) + builder.install_gradient_constraint() + assert "points" in seen + support = builder.interpolator.support + if sample_rule is None: + expected_local = support.barycentre + else: + from loop_interpolation._quadrature import tet_quadrature_points + + expected_local, _, _ = tet_quadrature_points(support.nodes, support.elements, sample_rule) + expected_world = model.bounding_box.reproject(expected_local) + np.testing.assert_allclose(seen["points"], expected_world, rtol=0, atol=1e-10) + world = expected_world + expected = r3_gradient(world, anchor) + expected /= np.linalg.norm(expected, axis=1)[:, None] + expected = model.bounding_box.project_vectors(expected) + np.testing.assert_allclose(seen["vectors"], expected, rtol=0, atol=1e-10) + + +@pytest.mark.parametrize("engine", ["FDI", "PLI"]) +@pytest.mark.parametrize("lower", [True, False]) +def test_inequality_feature_bounds_world(frame_case, engine, lower): + anchor, _, model = frame_case + builder = _builder(engine, model) + builder.interpolator.add_inequality_feature(_field(anchor), lower=lower) + bounds = builder.interpolator.ineq_constraints["inequality_feature"]["bounds"] + world = model.bounding_box.reproject(builder.interpolator.support.nodes) + values = r3_value(world, anchor) + column = 0 if lower else 1 + np.testing.assert_allclose(bounds[:, column], values, rtol=0, atol=tp(200, 200)) + + +@pytest.mark.parametrize("engine", ["FDI", "PLI"]) +def test_restricted_reference_exclusion_counts(frame_case, engine): + anchor, _, model = frame_case + builder = _builder(engine, model) + nodes = builder.interpolator.support.nodes + world = model.bounding_box.reproject(nodes) + lower, upper = anchor + (2.1, 2.1, 2.1), anchor + (5.9, 7.9, 9.9) + ref = _field(anchor) + ref.name = "restricted" + ref.support_box = (lower, upper) + selected = np.all((world >= lower) & (world <= upper), axis=1) + builder.add_equality_constraints(ref, lambda p: np.ones(len(p), dtype=bool)) + builder.install_equality_constraints() + equality = builder.interpolator.equal_constraints + equality_rows = np.concatenate([entry["col"] for entry in equality.values()]) + np.testing.assert_array_equal(np.sort(equality_rows), np.flatnonzero(selected)) + assert len(equality_rows) == int(selected.sum()) + assert "restricted_equality" in equality + + builder.add_orthogonal_feature(ref, B=0.0, sample_rule=None) + builder.install_gradient_constraint() + sample_local = builder.interpolator.support.barycentre + sample_world = model.bounding_box.reproject(sample_local) + valid_samples = np.all((sample_world >= lower) & (sample_world <= upper), axis=1) + orthogonal = builder.interpolator.constraints["restricted_orthogonal"] + assert orthogonal["matrix"].shape[0] == int(valid_samples.sum()) + excluded = builder.excluded_rows + assert excluded["restricted_equality"].get("invalid_reference_value", 0) == int( + (~selected).sum() + ) + assert excluded["restricted_orthogonal"].get("invalid_reference_direction", 0) == int( + (~valid_samples).sum() + ) + + +@pytest.mark.parametrize("engine", ["FDI", "PLI"]) +def test_support_expansion_contains_restored_points(frame_case, engine): + _, _, model = frame_case + builder = _builder(engine, model) + + class Restorer: + def __init__(self): + self.received = None + + def apply_to_points(self, points): + self.received = points.copy() + return points + np.array([20.0, -3.0, 2.0]) + + restorer = Restorer() + builder.faults.append(restorer) + model.bounding_box.nsteps = [2, 2, 2] + expected = model.bounding_box.regular_grid(local=False) + builder.check_interpolation_geometry({}, buffer=0) + np.testing.assert_allclose(restorer.received, expected, rtol=0, atol=tp(50, 50)) + restored = expected + [20.0, -3.0, 2.0] + local = model.bounding_box.project(restored) + support = builder.interpolator.support + assert np.all(local >= support.origin - 1e-10) + assert np.all(local <= support.maximum + 1e-10) + + +@pytest.mark.parametrize("engine", ["FDI", "PLI"]) +def test_set_interpolation_geometry_encloses_rotated_box(frame_case, engine): + anchor, rotation, model = frame_case + builder = _builder(engine, model) + corners = np.array(list(itertools.product(*zip(anchor, anchor + (8, 10, 12))))) + local = model.bounding_box.project(corners) + builder.set_interpolation_geometry(anchor, anchor + (8, 10, 12), rotation=rotation) + support = builder.interpolator.support + assert np.all(local >= support.origin - 1e-10) + assert np.all(local <= support.maximum + 1e-10)