fix: add caching of reused values

Author: Lachlan Grose

LoopStructural/modelling/features/_lambda_geological_feature.py

@@ -66,8 +66,25 @@ def evaluate_value(self, pos: np.ndarray, ignore_regions=False) -> np.ndarray:
         v = np.zeros((pos.shape[0]))
         v[:] = np.nan
 
-        mask = self._calculate_mask(pos, ignore_regions=ignore_regions)
-        pos = self._apply_faults(pos)
+        # Precompute each fault's scalar value (gx = fault.__getitem__(0).evaluate_value)
+        # once and reuse for both the region mask check and fault application.
+        # FaultSegment.evaluate_value(pos) == self.__getitem__(0).evaluate_value(pos) == gx.
+        # apply_to_points also needs gx to determine the displacement mask — passing it
+        # avoids the duplicate evaluation on the np.copy of pos.
+        precomputed_gx = {id(f): f.evaluate_value(pos) for f in self.faults}
+
+        # Region mask: pass precomputed gx so PositiveRegion/NegativeRegion skip re-evaluation
+        mask = np.ones(pos.shape[0], dtype=bool)
+        if not ignore_regions:
+            for r in self.regions:
+                pre = precomputed_gx.get(id(getattr(r, 'feature', None)))
+                mask = np.logical_and(mask, r(pos, precomputed_val=pre))
+
+        # Apply faults: pass precomputed gx so apply_to_points skips one evaluate_value call
+        if self.faults_enabled:
+            for f in self.faults:
+                pos = f.apply_to_points(pos, precomputed_gx=precomputed_gx.get(id(f)))
+
         if self.function is not None:
             v[mask] = self.function(pos[mask,:])
         return v

LoopStructural/modelling/features/fault/_fault_segment.py

@@ -311,13 +311,14 @@ def evaluate_displacement(self, points):
             d[mask] = self.faultfunction(gx[mask] + self.fault_offset, gy[mask], gz[mask])
         return d * self.displacement
 
-    def apply_to_points(self, points, reverse=False):
+    def apply_to_points(self, points, reverse=False, precomputed_gx=None):
         """
         Unfault the array of points
 
         Parameters
         ----------
         points - numpy array Nx3
+        precomputed_gx - optional pre-evaluated gx values (same points, avoids duplicate evaluation)
 
         Returns
         -------
@@ -328,10 +329,12 @@ def apply_to_points(self, points, reverse=False):
         newp = np.copy(points).astype(float)
         # evaluate fault function for all points
         # then define mask for only points affected by fault
-        gx = None
-        gy = None
-        gz = None
-        if use_threads:
+        # gx may be supplied by caller to avoid re-evaluation (precomputed from region check)
+        if precomputed_gx is not None:
+            gx = precomputed_gx
+            gy = self.__getitem__(1).evaluate_value(newp)
+            gz = self.__getitem__(2).evaluate_value(newp)
+        elif use_threads:
             with ThreadPoolExecutor(max_workers=8) as executor:
                 # all of these operations should be
                 # independent so just run as different threads
@@ -361,27 +364,32 @@ def apply_to_points(self, points, reverse=False):
             d *= -1.0
         # calculate the fault frame for the evaluation points
         for _i in range(steps):
-            gx = None
-            gy = None
-            gz = None
             g = None
-            if use_threads:
+            if _i == 0:
+                # Reuse gx/gy/gz from the initial full-array evaluation above — newp[mask] hasn't
+                # moved yet on the first iteration, so values are identical.
+                gx_m = gx[mask]
+                gy_m = gy[mask]
+                gz_m = gz[mask]
+                g = self.__getitem__(1).evaluate_gradient(newp[mask, :], ignore_regions=True)
+            elif use_threads:
                 with ThreadPoolExecutor(max_workers=8) as executor:
                     # all of these operations should be
                     # independent so just run as different threads
                     gx_future = executor.submit(self.__getitem__(0).evaluate_value, newp[mask, :])
                     g_future = executor.submit(self.__getitem__(1).evaluate_gradient, newp[mask, :])
                     gy_future = executor.submit(self.__getitem__(1).evaluate_value, newp[mask, :])
                     gz_future = executor.submit(self.__getitem__(2).evaluate_value, newp[mask, :])
-                    gx = gx_future.result()
+                    gx_m = gx_future.result()
                     g = g_future.result()
-                    gy = gy_future.result()
-                    gz = gz_future.result()
+                    gy_m = gy_future.result()
+                    gz_m = gz_future.result()
             else:
-                gx = self.__getitem__(0).evaluate_value(newp[mask, :])
-                gy = self.__getitem__(1).evaluate_value(newp[mask, :])
-                gz = self.__getitem__(2).evaluate_value(newp[mask, :])
+                gx_m = self.__getitem__(0).evaluate_value(newp[mask, :])
+                gy_m = self.__getitem__(1).evaluate_value(newp[mask, :])
+                gz_m = self.__getitem__(2).evaluate_value(newp[mask, :])
                 g = self.__getitem__(1).evaluate_gradient(newp[mask, :], ignore_regions=True)
+            gx, gy, gz = gx_m, gy_m, gz_m  # alias for block below
             # # get the fault frame val/grad for the points
             # determine displacement magnitude, for constant displacement
             # hanging wall should be > 0

LoopStructural/utils/regions.py

@@ -45,23 +45,22 @@ def __init__(self, feature, vector=None, point=None):
         self.name = 'PositiveRegion'
         self.parent = feature
 
-    def _calculate_value_and_distance(self, xyz)-> Tuple[np.ndarray, np.ndarray]:
-        val = self.feature.evaluate_value(xyz)
-        # find a point on/near 0 isosurface
+    def _calculate_value_and_distance(self, xyz, precomputed_val=None)-> Tuple[np.ndarray, np.ndarray]:
+        val = precomputed_val if precomputed_val is not None else self.feature.evaluate_value(xyz)
+        # find a point on/near 0 isosurface — compute once and cache on self
         if self.point is None:
             mask = np.zeros(xyz.shape[0], dtype="bool")
             mask[:] = val < 0
             if np.sum(mask) == 0:
                 raise ValueError("Cannot find point on surface")
-            centre = xyz[mask, :][0, :]
-        else:
-            centre = self.point
+            self.point = xyz[mask, :][0, :]  # cache: characterises the fault, not the query points
+        centre = self.point
         if self.vector is None:
             average_gradient = self.feature.evaluate_gradient(np.array([centre]))[0]
             average_gradient[2] = 0
             average_gradient /= np.linalg.norm(average_gradient)
-        else:
-            average_gradient = self.vector
+            self.vector = average_gradient  # cache: fault geometry doesn't change between calls
+        average_gradient = self.vector
 
         distance = np.einsum(
             "ij,j->i", centre[None, :] - xyz, average_gradient.reshape(-1, 3)[0, :]
@@ -80,8 +79,8 @@ def __init__(self, feature, vector=None, point=None):
         self.name = 'PositiveRegion'
         self.parent = feature
 
-    def __call__(self, xyz) -> np.ndarray:
-        val, distance = self._calculate_value_and_distance(xyz)
+    def __call__(self, xyz, precomputed_val=None) -> np.ndarray:
+        val, distance = self._calculate_value_and_distance(xyz, precomputed_val=precomputed_val)
         return np.logical_or(
             np.logical_and(~np.isnan(val), val > 0),
             np.logical_and(np.isnan(val), distance > 0),
@@ -99,8 +98,8 @@ def __init__(self, feature, vector=None, point=None):
         self.name = 'NegativeRegion'
         self.parent = feature
 
-    def __call__(self, xyz) -> np.ndarray:
-        val, distance = self._calculate_value_and_distance(xyz)
+    def __call__(self, xyz, precomputed_val=None) -> np.ndarray:
+        val, distance = self._calculate_value_and_distance(xyz, precomputed_val=precomputed_val)
         return np.logical_or(
             np.logical_and(~np.isnan(val), val < 0),
             np.logical_and(np.isnan(val), distance < 0),