fix: point transformations, fault-frame validation and backend matrix, stage S6 (spec 005)

COMPAT.md

@@ -56,11 +56,20 @@ unchanged, but the new path is preferred going forward.
 | `GeologicalFeature.get_data()` (used by `Loop2DView.add_data` and `ModelExporter.save`) returned locations moved by the local origin and rotated a second time; the intrusion classification evaluated the model at such points | Both use the stored world locations. With an identity frame nothing changes (spec 005, stage S4) | 2026-11-05 |
 | `GeologicalFeature.evaluate_gradient_misfit()` changed the stored constraint arrays in place, compared a local vector with a world vector and read one column. Its value was a dot product with an unnormalised model gradient | The name and the public shape stay. The method changes no stored array. It gives the cosine between the world data direction and the model gradient for each row, so the value differs from the old value. Use `evaluate_gradient_cosines()` for a 1-D array. Two families of different length give one row of shape `(1, N1 + N2)`; before, this raised `ValueError` (spec 005, stage S5) | 2026-11-05 |
 | `evaluate_gradient_misfit()` and `evaluate_value_misfit()` rebuilt a stale model | Both methods write one `UserWarning` and return NaN. They do not rebuild. Build the model first. The new public method `GeologicalFeature.verification_report()` gives the constraint diagnostics (residuals and exclusions) and does not build the model. A PLI model ignores gradient data, so its report shows an `unattributed` record for the gradient rows (spec 005, stage S5) | 2026-11-05 |
-| `evaluate_value_misfit()` divided by the field range, with a division by zero for a zero range | It keeps the normalisation by the field range. A zero or not finite range gives the absolute difference. `evaluate_value_residuals()` gives the unnormalised 1-D result (spec 005, stage S5) | 2026-11-05 |
+| `evaluate_value_misfit()` divided by the field range, with a division by zero for a zero range | It keeps the normalisation by the field range. A zero or not finite range gives the absolute difference. `evaluate_value_residuals()` gives the unnormalised 1-D result. With an identity frame and no resize, the value misfit is the same as before (spec 005, stage S5) | 2026-11-05 |
 | Assignment of `nsteps` or `step_vector` on a 3D structured grid, and `set_nelements` | The extent stays the same. Assignment gives a new step, not a product rule. `nsteps` must be whole and at least 2. `set_nelements` keeps `maximum`. Grid steps and solved values can change slightly (spec 005, stage S5) | 2026-11-05 |
 | Assignment of `origin`, `nsteps` or `step_vector` on a `RectilinearGrid` | It raises `NotImplementedError`. The node arrays are fixed (spec 005, stage S5) | 2026-11-05 |
 | The 2D structured grid: `nsteps_cells`, `n_cell_x` and `n_cell_y` were separate values, and `maximum` could be assigned | The grid derives `nsteps_cells`, `n_cell_x` and `n_cell_y` from `nsteps`. `maximum` is read-only (assignment raises `AttributeError`) and keeps its legacy value (spec 005, stage S5) | 2026-11-05 |
-| `ConstraintDiagnosticsReport` had no residual or exclusion fields | It has three new optional fields: `physical_residuals`, `assembled_residuals` and `exclusions`. With an identity frame and no resize, the value misfit is the same as before (spec 005, stage S5) | 2026-11-05 |
+| `ConstraintDiagnosticsReport` had no residual or exclusion fields | It has three new optional fields: `physical_residuals`, `assembled_residuals` and `exclusions` (spec 005, stage S5) | 2026-11-05 |
+| `LoopStructural.utils.EuclideanTransformation` was a separate class in `LoopStructural`. Its `inverse_transform` failed for an `(N, d)` array (`ValueError` or too few rows), its methods truncated an integer input, and a 1-D input raised `IndexError` | `EuclideanTransformation` is one class. `loop_common.math.EuclideanTransformation` is new, and `LoopStructural.utils.EuclideanTransformation` is the same object. `inverse_transform` returns the original points for an `(N, d)` array. The methods return a float copy. A 1-D input raises `ValueError` (spec 005, stage S6) | 2026-11-06 |
+| `EuclideanTransformation` with `dimensions=3`: z was reversed, `fit` failed, and the fitted angle had the wrong sign for a main axis with a negative slope (a line with a direction between -45 and 0 degrees). `fit` with `fit_rotation=True` logged an error and returned `None` when scikit-learn was missing | For `dimensions=3` z is no longer reversed (`z' = z - tz`; before, `-(z - tz)`). `fit` works for `dimensions=3`. The fitted angle has the right sign for a main axis with a negative slope. A 3D user must check every result of the old class. `fit` with `fit_rotation=True` raises `ImportError` for 2 or more rows when scikit-learn is missing. With fewer than 2 rows the angle is 0 and scikit-learn is not used. `fit` with `fit_rotation=False` needs no scikit-learn (spec 005, stage S6) | 2026-11-06 |
+| `EuclideanTransformation.fit` logged through the `LoopStructural` logger | `fit` logs through the `loop_common` logger (`loop_common.logging.get_logger`) (spec 005, stage S6) | 2026-11-06 |
+| A pickle that holds an `EuclideanTransformation` | A pickle made before S6 still loads. The old module path stays importable and exposes the class. A pickle made after S6 stores `loop_common.math._transformation`. A release from before S6 cannot load it (spec 005, stage S6) | 2026-11-06 |
+| `evaluate_gradient_misfit`, `evaluate_value_misfit`, `evaluate_gradient_cosines`, `evaluate_value_residuals` and `verification_report` of a faulted feature compared restored rows with the interpolator after the faults acted a second time. For a feature with a region mask, a row outside the region gave NaN | These methods compare the stored restored rows with the interpolator at the same positions. `_value_residuals_fresh` and `_gradient_cosines_fresh` now evaluate the interpolator directly. This reverses S5 decision 8 for these two methods. Numbers for a feature with no faults and no regions do not change. A row outside a region mask now gives a number (spec 005, stage S6) | 2026-11-06 |
+| Known issue, not fixed in S6: the public `evaluate_gradient` of a faulted feature | It keeps its tetrahedron path and the shape error that S3 noted (spec 005, stage S6) | 2026-11-06 |
+| Known issue, not fixed in S6: the gradient of a model outside the support | Outside the support the FDI gradient is 0 and the PLI gradient is NaN. S6 test T095 pins this behavior (spec 005, stage S6) | 2026-11-06 |
+| Known issue, not fixed in S6: `P2Interpolator.add_gradient_constraints` for a P2 foliation with gradient rows | It raises `ValueError`. This is a library defect. S6 does not change it (spec 005, stage S6) | 2026-11-06 |
+| Known issue, not fixed in S6: `set_nelements` | It passes a cell count to the node setter. For example, `nelements=1000` gives 729 cells. This is old behavior, not a regression. S6 does not change it (spec 005, stage S6, S5 finding 3) | 2026-11-06 |
 
 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.
 

LoopStructural/modelling/features/_geological_feature.py

@@ -334,6 +334,9 @@ def _gradient_cosines_fresh(self) -> list:
         """Compute the gradient and normal cosines without a stale check.
 
         The method copies the stored arrays. It does not change them.
+        The stored positions are already restored by the faults of the feature.
+        The method evaluates the interpolator, not the feature, so no fault
+        acts a second time.
 
         Returns
         -------
@@ -355,13 +358,17 @@ def _gradient_cosines_fresh(self) -> list:
             if bounding_box is not None:
                 positions = bounding_box.reproject(positions)
                 data = self.interpolator.local_to_world_vectors(data)
-            model = self.evaluate_gradient(positions)
+            model = self.interpolator.evaluate_gradient(positions)
             cosines.append(direction_cosines(model, data))
         return cosines
 
     def _value_residuals_fresh(self) -> np.ndarray:
         """Compute the value residuals without a stale check.
 
+        The stored positions are already restored by the faults of the feature.
+        The method evaluates the interpolator, not the feature, so no fault
+        acts a second time.
+
         Returns
         -------
         numpy.ndarray
@@ -371,7 +378,7 @@ def _value_residuals_fresh(self) -> np.ndarray:
         positions = rows[:, :3]
         if self.interpolator.bounding_box is not None:
             positions = self.interpolator.bounding_box.reproject(positions)
-        return np.abs(rows[:, 3] - self.evaluate_value(positions))
+        return np.abs(rows[:, 3] - self.interpolator.evaluate_value(positions))
 
     def _physical_residual_summaries(self) -> dict[str, ResidualSummary]:
         """Summarise the physical residuals of the value, gradient and normal data.

LoopStructural/utils/_transformation.py

@@ -1,176 +1,8 @@
-import numpy as np
+"""Old import path of ``EuclideanTransformation``.
 
-from .logging import getLogger
+The class lives in ``loop_common.math``. This module keeps the path
+``LoopStructural.utils._transformation`` working.
+"""
+from loop_common.math._transformation import EuclideanTransformation
 
-logger = getLogger(__name__)
-
-
-class EuclideanTransformation:
-    def __init__(
-        self,
-        dimensions: int = 2,
-        angle: float = 0,
-        translation: np.ndarray = None,
-        fit_rotation: bool = True,
-    ):
-        """Transforms points into a new coordinate
-        system where the main eigenvector is aligned with x
-
-        Parameters
-        ----------
-        dimensions : int, optional
-            Do transformation in map view or on 3d volume, by default 2
-        angle : float, optional
-            Angle to rotate the points by, by default 0
-        translation : np.ndarray, default zeros
-            Translation to apply to the points, by default
-        """
-        if translation is None:
-            translation = np.zeros(3)
-        self.translation = translation[:dimensions]
-        self.dimensions = dimensions
-        self.angle = angle
-        self.fit_rotation = fit_rotation
-
-    def fit(self, points: np.ndarray):
-        """Fit the transformation to a point cloud
-        This function will find the main eigenvector of the point cloud
-        and rotate the point cloud so that this is aligned with x
-
-
-        Parameters
-        ----------
-        points : np.ndarray
-            xyz points as as numpy array
-        """
-        try:
-            from sklearn import decomposition
-        except ImportError:
-            logger.error('scikit-learn is required for this function')
-            return
-        points = np.array(points)
-        if points.shape[1] < self.dimensions:
-            raise ValueError(f"Points must have at least {self.dimensions} dimensions")
-        # standardise the points so that centre is 0
-        # self.translation = np.zeros(3)
-        self.translation = np.mean(points[:, : self.dimensions], axis=0)
-        # find main eigenvector and and calculate the angle of this with x
-        if self.fit_rotation:
-            pca = decomposition.PCA(n_components=self.dimensions).fit(
-                points[:, : self.dimensions] - self.translation[None, : self.dimensions]
-            )
-            coeffs = pca.components_
-            self.angle = -np.arccos(np.dot(coeffs[0, :], [1, 0]))
-        else:
-            self.angle = 0
-        return self
-
-    @property
-    def rotation(self):
-        return self._rotation(self.angle)
-
-    @property
-    def inverse_rotation(self):
-        return self._rotation(-self.angle)
-
-    def _rotation(self, angle):
-        return np.array(
-            [
-                [np.cos(angle), -np.sin(angle), 0],
-                [np.sin(angle), np.cos(angle), 0],
-                [0, 0, -1],
-            ]
-        )
-
-    def fit_transform(self, points: np.ndarray) -> np.ndarray:
-        """Fit the transformation and transform the points"""
-
-        self.fit(points)
-        return self.transform(points)
-
-    def transform(self, points: np.ndarray) -> np.ndarray:
-        """Transform points using the transformation and rotation
-
-        Parameters
-        ----------
-        points : np.ndarray
-            xyz points as as numpy array
-
-        Returns
-        -------
-        np.ndarray
-            xyz points in the transformed coordinate system
-        """
-        points = np.array(points)
-        if points.shape[1] < self.dimensions:
-            raise ValueError(f"Points must have at least {self.dimensions} dimensions")
-        centred = points[:, : self.dimensions] - self.translation[None, :]
-        rotated = np.einsum(
-            'ik,jk->ij',
-            centred,
-            self.rotation[: self.dimensions, : self.dimensions],
-        )
-        transformed_points = np.copy(points)
-        transformed_points[:, : self.dimensions] = rotated
-        return transformed_points
-
-    def inverse_transform(self, points: np.ndarray) -> np.ndarray:
-        """
-        Transform points back to the original coordinate system
-
-        Parameters
-        ----------
-        points : np.ndarray
-            xyz points as as numpy array
-
-        Returns
-        -------
-        np.ndarray
-            xyz points in the original coordinate system
-        """
-        inversed = (
-            np.einsum(
-                'ik,jk->ij',
-                points[: self.dimensions],
-                self.inverse_rotation[: self.dimensions, : self.dimensions],
-            )
-            + self.translation
-        )
-        inversed = (
-            np.vstack([inversed, points[self.dimensions :]])
-            if points.shape[1] > self.dimensions
-            else inversed
-        )
-        return inversed
-
-    def __call__(self, points: np.ndarray) -> np.ndarray:
-        """
-        Transform points into the transformed space
-
-        Parameters
-        ----------
-        points : np.ndarray
-            xyz points as as numpy array
-
-        Returns
-        -------
-        np.ndarray
-            xyz points in the transformed coordinate system
-        """
-
-        return self.transform(points)
-
-    def _repr_html_(self):
-        """
-        Provides an HTML representation of the TransRotator.
-        """
-        html_str = f"""
-        <div class="collapsible">
-          <button class="collapsible-button">{self.__class__.__name__}</button>
-          <div class="content">
-            <p>Translation: {self.translation}</p>
-            <p>Rotation Angle: {self.angle} degrees</p>
-          </div>
-        </div>
-        """
-        return html_str
+__all__ = ["EuclideanTransformation"]

ROADMAP.md

@@ -152,6 +152,26 @@ A stale model gives a `UserWarning` and NaN. The model is not rebuilt.
 The new methods `evaluate_gradient_cosines` and `evaluate_value_residuals`
 give the 1-D result. The notice is in `COMPAT.md`.
 
+### Residuals of faulted features (spec 005, S6)
+Decision date: 2026-11-06 (spec 005, S6).
+
+The residual diagnostics of a faulted feature evaluate the interpolator at
+the restored positions. S5 decision 8 changes for the two private methods
+`_value_residuals_fresh` and `_gradient_cosines_fresh`. The notice is in
+`COMPAT.md`.
+
+### EuclideanTransformation (spec 005, S6)
+Decision date: 2026-11-06 (spec 005, S6).
+
+This is a behaviour change. The class had four defects: the 3D `fit`, the
+sign of the fitted angle, the reversed z axis, and the silent need for
+scikit-learn. The physical effect was a mirrored z axis and a wrong rotation
+of the map view for some lines.
+
+One class now lives in `loop_common`. `LoopStructural.utils` re-exports
+it. This reverses item 2c-5 of the stage sequence, which was closed as
+deferred. The notice is in `COMPAT.md`.
+
 ## Stage sequence
 
 - [x] **Stage 0 — Planning infra.** This file, release/versioning/compat
@@ -431,6 +451,8 @@ give the 1-D result. The notice is in `COMPAT.md`.
     (`translation: np.ndarray = np.zeros(3)`) as part of the merge.
     Closed as deferred: local class remains the runtime source for now;
     mutable-default regression was already eliminated in LoopStructural.
+    Spec 005 S6 reverses this closure: the one class in `loop_common` is now
+    the runtime source (see the S6 decision in `## Decisions`).
   - [x] **2c-6.** Resolve `loop_common`'s reverse dependency on
     `LoopStructural.export.*`: move `LoopStructural/export/geoh5.py`,
     `gocad.py`, `omf_wrapper.py`, `exporters.py` into `loop_common/io/`
@@ -894,3 +916,4 @@ give the 1-D result. The notice is in `COMPAT.md`.
 - **2026-11-05:** Spec 005 stage S3 implemented: saved interpolator frame, legacy-load rule, frame-keeping copies, stale-frame check. Gate result: root pytest 1978 passed, 5 skipped.
 - **2026-11-05:** Spec 005 stage S4 implemented: display, export and intrusion classification use stored world locations; rescale and scale deprecated.
 - **2026-11-05:** Spec 005 stage S5 implemented: resize paths keep the extent, cell-centre guard tests, non-mutating world-correct misfit diagnostics, residual and exclusion report.
+- **2026-11-06:** Spec 005 complete: stages S1 to S6 implemented; the auxiliary transformation, the R2 displacement validation, the backend matrix and the golden review are done. Gate result: root pytest 2781 passed, 30 skipped (every skip names its dependency in `specs/005-coordinate-frame-correctness/verification.md`).

packages/loop_common/src/loop_common/geometry/_structured_grid_3d.py

@@ -1,4 +1,4 @@
-"""Pure 3D regular grid geometry: origin/nsteps/step_vector indexing."""
+"""Pure 3D regular grid geometry: origin/nsteps/step_vector indexing."""
 
 from typing import Optional
 

packages/loop_common/src/loop_common/math/__init__.py

@@ -1,5 +1,6 @@
 import numpy as np
 
 from ._maths import *
+from ._transformation import EuclideanTransformation
 
 rng = np.random.default_rng()