fix: Updating intrusions code

- Added a new method `_validate_intrusion_inputs` in `GeologicalModel` to validate inputs for intrusions, ensuring necessary data is present before processing.
- Updated `_build_intrusion` to call the new validation method, improving error handling for missing data.
- Refactored `IntrusionBuilder.create_geometry_using_geometric_scaling` to clarify that geometric scaling is not currently implemented, raising a `NotImplementedError` immediately.
- Simplified threshold handling in `IntrusionFeature` by removing redundant checks for marginal faults.
- Removed the unused `intrusion_support_functions.py` file to clean up the codebase.
- Updated tests in `test_intrusions.py` to cover new validation logic, ensuring clear error messages for missing data and parameters.
- Added regression tests for previously silent errors related to weight handling and geometric scaling.

(cherry picked from commit 8cb80c1)

Author: lachlangrose

INTRUSIONS.md

@@ -1402,6 +1402,50 @@ def create_and_add_intrusion(
             **kwargs,
         )
 
+    def _validate_intrusion_inputs(
+        self,
+        intrusion_name,
+        intrusion_frame_name,
+        intrusion_data,
+        intrusion_frame_data,
+        intrusion_frame_parameters,
+    ):
+        """Fail fast, at the `create_and_add_intrusion` boundary, with a clear
+        message naming the missing piece -- instead of a bare `KeyError`
+        several calls deep inside `IntrusionFrameBuilder`/`IntrusionBuilder`
+        once building has already started. See ``INTRUSIONS.md`` finding 4.
+        """
+        if intrusion_data.empty:
+            raise ValueError(
+                f"No data found for intrusion '{intrusion_name}': check that "
+                "model.data contains rows with feature_name == "
+                f"'{intrusion_name}'"
+            )
+        if intrusion_frame_data.empty:
+            raise ValueError(
+                f"No data found for intrusion frame '{intrusion_frame_name}': "
+                "check that model.data contains rows with feature_name == "
+                f"'{intrusion_frame_name}'"
+            )
+        required_columns = ["intrusion_contact_type", "intrusion_side"]
+        missing_columns = [c for c in required_columns if c not in intrusion_data.columns]
+        if missing_columns:
+            raise ValueError(
+                f"Intrusion data for '{intrusion_name}' is missing required "
+                f"column(s) {missing_columns}: 'intrusion_contact_type' marks "
+                "each point as 'roof'/'top' or 'floor'/'base', and "
+                "'intrusion_side' (boolean) marks points used to constrain "
+                "the lateral extent"
+            )
+        contact_anisotropies = intrusion_frame_parameters.get("contact_anisotropies")
+        if not contact_anisotropies:
+            raise ValueError(
+                "intrusion_frame_parameters['contact_anisotropies'] is "
+                "required: provide a non-empty list of series-type features "
+                "to use as the inflation-gradient proxy for the intrusion "
+                "frame's coordinate 0"
+            )
+
     def _build_intrusion(
         self,
         intrusion_name,
@@ -1458,6 +1502,14 @@ def _build_intrusion(
         intrusion_data = self.data[self.data["feature_name"] == intrusion_name].copy()
         intrusion_frame_data = self.data[self.data["feature_name"] == intrusion_frame_name].copy()
 
+        self._validate_intrusion_inputs(
+            intrusion_name,
+            intrusion_frame_name,
+            intrusion_data,
+            intrusion_frame_data,
+            intrusion_frame_parameters,
+        )
+
         # -- get variables for intrusion frame interpolation
         gxxgz = kwargs.get("gxxgz", 0)
         gxxgy = kwargs.get("gxxgy", 0)
@@ -1496,7 +1548,7 @@ def _build_intrusion(
             nelements=nelements,
             w2=weights[0],
             w1=weights[1],
-            gxygz=weights[2],
+            gyxgz=weights[2],
         )
 
         intrusion_frame = intrusion_frame_builder.frame

LoopStructural/modelling/core/geological_model.py

@@ -3,7 +3,6 @@
 
 from ...utils import getLogger, rng
 from ..features.builders import BaseBuilder
-from .geometric_scaling_functions import *
 from .intrusion_feature import IntrusionFeature
 
 logger = getLogger(__name__)
@@ -113,45 +112,25 @@ def set_data_for_extent_calculation(self, intrusion_data: pd.DataFrame):
     def create_geometry_using_geometric_scaling(
         self, geometric_scaling_parameters, reference_contact_data
     ):
-
-        geometric_scaling_parameters.get("intrusion_type", None)
-        intrusion_length = geometric_scaling_parameters.get("intrusion_length", None)
-        geometric_scaling_parameters.get("inflation_vector", np.array([[0, 0, 1]]))
-        thickness = geometric_scaling_parameters.get("thickness", None)
-
-        if (
-            self.intrusion_frame.builder.intrusion_network_contact == "floor"
-            or self.intrusion_frame.builder.intrusion_network_contact == "base"
-        ):
-            geometric_scaling_parameters.get("inflation_vector", np.array([[0, 0, 1]]))
-        else:
-            geometric_scaling_parameters.get("inflation_vector", np.array([[0, 0, -1]]))
-
-        if intrusion_length is None and thickness is None:
-            raise ValueError(
-                f"No {self.intrusion_frame.builder.intrusion_other_contact} data. Add intrusion_type and intrusion_length (or thickness) to geometric_scaling_parameters dictionary"
-            )
-
-        else:  # -- create synthetic data to constrain interpolation using geometric scaling
-            estimated_thickness = thickness
-            if estimated_thickness is None:
-                raise NotImplementedError("Not implemented")
-                # estimated_thickness = thickness_from_geometric_scaling(
-                #     intrusion_length, intrusion_type
-                # )
-
-            logger.info(
-                f"Building tabular intrusion using geometric scaling parameters: estimated thicknes = {round(estimated_thickness)} meters"
-            )
-            raise NotImplementedError("Not implemented")
-            # (
-            #     other_contact_data_temp,
-            #     other_contact_data_xyz_temp,
-            # ) = contact_pts_using_geometric_scaling(
-            #     estimated_thickness, reference_contact_data, inflation_vector
-            # )
-
-            # return other_contact_data_temp
+        """Not currently implemented.
+
+        This is meant to synthesise the missing contact (roof or floor) from
+        an estimated thickness (either given directly or derived from
+        empirical length/thickness scaling laws, see
+        ``geometric_scaling_functions.thickness_from_geometric_scaling``) and
+        an inflation vector, via
+        ``geometric_scaling_functions.contact_pts_using_geometric_scaling``.
+        That wiring was never completed, so every call path here always
+        raised ``NotImplementedError`` regardless of what was passed in
+        (see ``INTRUSIONS.md`` finding 2). Raising immediately, rather than
+        after partially validating parameters, makes that unambiguous.
+        """
+        raise NotImplementedError(
+            "geometric_scaling_parameters is not currently supported: "
+            f"'{self.intrusion_frame.builder.intrusion_other_contact}' contact "
+            "has no data, and synthesising it from geometric scaling is not "
+            "implemented. Provide explicit data for both contacts instead."
+        )
 
     def prepare_data(self, geometric_scaling_parameters):
         """Prepare the data to compute distance thresholds along the frame coordinates.

LoopStructural/modelling/intrusions/intrusion_builder.py

@@ -265,13 +265,8 @@ def evaluate_value(self, pos):
             intrusion_coord1_pts
         )
 
-        if self.intrusion_frame.builder.marginal_faults is not None:
-            c2_minside_threshold = thresholds[0]  # np.zeros_like(intrusion_coord2_pts)
-            c2_maxside_threshold = thresholds[1]
-
-        else:
-            c2_minside_threshold = thresholds[0]
-            c2_maxside_threshold = thresholds[1]
+        c2_minside_threshold = thresholds[0]
+        c2_maxside_threshold = thresholds[1]
 
         thresholds, _residuals, _conceptual = self.interpolate_vertical_thresholds(
             intrusion_coord1_pts, intrusion_coord2_pts
@@ -332,81 +327,6 @@ def evaluate_value(self, pos):
 
         return intrusion_sf
 
-    def evaluate_value_test(self, points):
-        """
-        Computes a distance scalar field to the intrusion contact (isovalue = 0).
-
-        Parameters
-        ------------
-        points : numpy array (x,y,z),  points where the IntrusionFeature is evaluated.
-
-        Returns
-        ------------
-        intrusion_sf : numpy array, contains distance to intrusion contact
-
-        """
-        self.builder.up_to_date()
-
-        # compute coordinates values for each evaluated point
-        intrusion_coord0_pts = self.intrusion_frame[0].evaluate_value(points)
-        intrusion_coord1_pts = self.intrusion_frame[1].evaluate_value(points)
-        intrusion_coord2_pts = self.intrusion_frame[2].evaluate_value(points)
-
-        self.evaluated_points = [
-            points,
-            intrusion_coord0_pts,
-            intrusion_coord1_pts,
-            intrusion_coord2_pts,
-        ]
-
-        thresholds, _residuals, _conceptual = self.interpolate_lateral_thresholds(
-            intrusion_coord1_pts
-        )
-
-        if self.intrusion_frame.builder.marginal_faults is not None:
-            c2_minside_threshold = np.zeros_like(intrusion_coord2_pts)
-            c2_maxside_threshold = thresholds[1]
-
-        else:
-            c2_minside_threshold = thresholds[0]
-            c2_maxside_threshold = thresholds[1]
-
-        thresholds, _residuals, _conceptual = self.interpolate_vertical_thresholds(
-            intrusion_coord1_pts, intrusion_coord2_pts
-        )
-        c0_minside_threshold = thresholds[1]
-        c0_maxside_threshold = thresholds[0]
-
-        mid_point = c0_minside_threshold + ((c0_maxside_threshold - c0_minside_threshold) / 2)
-
-        mod_intrusion_coord0_pts = intrusion_coord0_pts - mid_point
-        mod_c0_minside_threshold = c0_minside_threshold - mid_point
-        mod_c0_maxside_threshold = c0_maxside_threshold + mid_point
-
-        a = (
-            (mod_intrusion_coord0_pts >= mid_point)
-            * (c2_minside_threshold < intrusion_coord2_pts)
-            * (intrusion_coord2_pts < c2_maxside_threshold)
-        )
-        b = (
-            (mod_intrusion_coord0_pts <= mid_point)
-            * (c2_minside_threshold < intrusion_coord2_pts)
-            * (intrusion_coord2_pts < c2_maxside_threshold)
-        )
-        c = (
-            (mod_intrusion_coord0_pts <= mid_point)
-            * (mod_intrusion_coord0_pts >= mod_c0_minside_threshold)
-            * (c2_minside_threshold < intrusion_coord2_pts)
-            * (intrusion_coord2_pts < c2_maxside_threshold)
-        )
-
-        intrusion_sf = mod_intrusion_coord0_pts
-        intrusion_sf[a] = mod_intrusion_coord0_pts[a] - mod_c0_maxside_threshold[a]
-        intrusion_sf[b] = abs(mod_c0_minside_threshold[b] + mod_intrusion_coord0_pts[b])
-        intrusion_sf[c] = mod_intrusion_coord0_pts[c] - mod_c0_minside_threshold[c]
-
-        return intrusion_sf
-
     def get_data(self, value_map: dict | None = None):
         pass
 

LoopStructural/modelling/intrusions/intrusion_feature.py

@@ -17,6 +17,12 @@
     logger.error('Scikitlearn cannot be imported')
     raise
 
+# Fixed (not derived from the shared `rng`) so that repeated builds of the
+# same intrusion produce the same contact/fault clustering: `loop_common`'s
+# shared `rng` is a fresh, unseeded `np.random.default_rng()` per process, so
+# routing this through it would make cluster labels vary run-to-run instead.
+_KMEANS_RANDOM_STATE = 0
+
 
 class IntrusionFrameBuilder(StructuralFrameBuilder):
     def __init__(
@@ -228,10 +234,9 @@ def add_contact_anisotropies(self, series_list: list | None = None, **kwargs):
 
             # -- use scalar field values to find different contacts
             series_i_vals_mod = series_i_vals.reshape(len(series_i_vals), 1)
-            # TODO create global loopstructural random state variable
-            contact_clustering = KMeans(n_clusters=n_contacts, random_state=0).fit(
-                series_i_vals_mod
-            )
+            contact_clustering = KMeans(
+                n_clusters=n_contacts, random_state=_KMEANS_RANDOM_STATE
+            ).fit(series_i_vals_mod)
 
             for j in range(n_contacts):
                 z = np.ma.masked_not_equal(contact_clustering.labels_, j)
@@ -358,7 +363,9 @@ def set_intrusion_steps_parameters(self):
                 )
                 series_values = series_from_name.evaluate_value(data_points_xyz)
                 series_values_mod = series_values.reshape(len(series_values), 1)
-                contact_clustering = KMeans(n_clusters=2, random_state=0).fit(series_values_mod)
+                contact_clustering = KMeans(
+                    n_clusters=2, random_state=_KMEANS_RANDOM_STATE
+                ).fit(series_values_mod)
 
                 # contact 0
                 z = np.ma.masked_not_equal(contact_clustering.labels_, 0)
@@ -491,7 +498,6 @@ def set_marginal_faults_parameters(self):
         for fault_i in self.marginal_faults:
             marginal_fault = self.marginal_faults[fault_i].get("structure")
             block = self.marginal_faults[fault_i].get("block")  # hanging wall or foot wall
-            self.marginal_faults[fault_i].get("emplacement_mechanism")
             series_name = self.marginal_faults[fault_i].get("series")
 
             series_values_temp = series_name.evaluate_value(intrusion_frame_c0_data_xyz)
@@ -828,7 +834,6 @@ def create_constraints_for_c0(self, **kwargs):
                 delta_contact = self.marginal_faults[fault_i].get("delta_c", 1)
                 marginal_fault = self.marginal_faults[fault_i].get("structure")
                 block = self.marginal_faults[fault_i].get("block")  # hanging wall or foot wall
-                self.marginal_faults[fault_i].get("emplacement_mechanism")
                 series_name = self.marginal_faults[fault_i].get("series")
 
                 fault_gridpoints_vals = marginal_fault[0].evaluate_value(grid_points)
@@ -976,7 +981,3 @@ def set_intrusion_frame_data(self, intrusion_frame_data):  # , intrusion_network
 
         self.add_data_from_data_frame(intrusion_frame_data_complete)
         self.update_geometry(intrusion_frame_data_complete[["X", "Y", "Z"]].to_numpy())
-
-    def update(self):
-        for i in range(3):
-            self.builders[i].update()